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ABSTRACT 


In  this  thesis,  we  study  the  artifacts  that  occur  when  a  scene  being  imaged  by  radar 
contains  moving  targets.  The  physics  of  interaction  between  radar  waves  and  moving 
targets  were  studied  to  develop  a  model  using  MATLAB  for  the  received  signal,  which 
does  not  make  use  of  the  start- stop  approximation.  The  effects  of  target  motion  on  the 
image  formation  process  were  studied  for  different  radar  configurations,  including 
multistatic  radars  and  Synthetic  Aperture  Radars.  The  key  limitation  of  this  model  is  its 
high  computational  resource  requirements  when  simulating  high  bandwidth  or  long 
pulses. 

It  was  observed  that  range  profiles  may  experience  distortion  due  to  the  received 
waveform’s  differences  from  the  matched  filter.  The  exact  outcome  is  waveform 
dependent;  generally,  both  main  lobe  broadening  and  range  errors  were  introduced  by 
target  motion.  This  leads  to  incorrect  object  localization  and  defocusing  on  the  image. 
For  SAR,  a  moving  target’s  physical  location  varies  throughout  the  imaging  process.  This 
means  that  standard  backprojection  fails  to  yield  a  focused  image  even  if  the  range  error 
due  to  the  Doppler  shift  has  been  corrected,  resulting  in  smearing.  This  is  similar  to  the 
“motion  blur”  experienced  in  optical  cameras  with  a  fast  object. 
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I.  INTRODUCTION 


A.  HISTORY  AND  PRINCIPLES  OF  RADAR 

The  physical  principles  behind  radar  operations  are  often  said  to  have  been  first 
demonstrated  by  Heinrich  Hertz  in  his  series  of  experiments  on  the  properties  of 
electromagnetic  waves  in  the  late  1 880— 1 890s.  In  April  1904,  a  German  Engineer 
Christian  Hiilsmeyer  patented  a  device  in  Germany  for  ship  collision  avoidance  called  a 
“Telemobiloskop,”  which  had  the  complete  functionality  of  radar  [1],  Thus,  the  idea 
behind  radar — which  is  an  acronym  for  “RAdio  Detection  And  Ranging” — has  been 
around  for  more  than  a  hundred  years.  The  basic  working  principles  of  radars  can  be 
simply  summarized  as  follows: 

•  A  radio  transmitter  sends  out  a  known  electromagnetic  wave  that  travels 
through  the  surrounding  transmission  media  (e.g.,  the  atmosphere) 

•  This  wave  reaches  and  interacts  with  some  objects.  Part  of  the  wave  may 
be  reflected/deflected  in  different  directions. 

•  A  radio  receiver  senses  the  portion  of  the  wave  that  emanates  back  from 
the  objects  and  compares  the  characteristics  (reception  time,  frequency, 
phase,  etc.)  of  the  received  wave  versus  the  transmitted  wave  to  extract 
information  (position,  speed,  etc.)  about  the  object  through  processing  of 
the  received  signal. 

The  prime  advantage  of  radar  over  other  forms  of  detection  utilizing  the  visible 
light  and  infra-red  portions  of  electromagnetic  spectrum  lies  in  the  fact  that  the 
microwaves  that  radars  utilize  are  relatively  unaffected  by  weather  or  atmospheric 
absorption. 

Even  amongst  radars,  there  are  many  different  types.  The  electromagnetic  wave 
sent  out  by  a  radar  may  either  be  pulses  of  energy,  leading  to  what  is  called  “pulsed  radar,” 
or  they  may  be  a  continuous  transmission  of  energy,  leading  to  what  are  called 
“continuous  wave”  or  “CW  radars.”  Radars  where  the  transmitters  and  receivers  are  at 
different  locations  are  termed  “multi-static;”  and  they  are  termed  “mono-static”  if  these 
are  co-located.  The  work  done  here  addresses  pulsed  radars  in  both  mono-static  and 
multi-static  configurations. 
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Figure  1  illustrates  the  key  principles  and  phenomena  that  govern  the  operation  of 
a  typical  radar. 
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Figure  1.  Key  principles  and  phenomena  for  operation  of  typical  pulse  radar 

(After  [2]) 


B.  RADAR  IMAGING 

The  concept  of  radar  imaging  represents  an  extension  of  the  conventional  use  of 
radar  and  was  an  idea  developed  from  the  1950’s  onwards  following  the  realization  of  its 
possibilities  by  Carl  Wiley  of  Goodyear  Aircraft  Corporation.  Radar  imaging  basically 
seeks  to  form  a  spatial  image  of  the  radar  reflectivity  of  an  observed  scene  or  target  using 
an  analysis  of  scattered  wave  returns.  While  the  intent,  like  optical  images,  is  to  draw 
information  from  the  imaged  scene,  radar  imagery  typically  looks  different  from  the 
optical  equivalent.  Figure  2  shows  an  example  of  a  radar  image  versus  an  optical  image. 
At  radar  wavelengths,  the  backscattering  of  waves  depends  on  different  factors  than  those 
of  optical  wavelengths.  In  radar  imagery,  the  surface  roughness  determines  how  much  of 
the  wave  is  backscattered  to  the  receiver;  in  general,  a  rougher  surface  backscatters  more 
of  the  wave  to  the  receiver.  As  the  backscattered  waves  are  coherently  processed  in 
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radars,  objects  where  the  waves  coherently  sum  up  (e.g.,  comer  reflectors)  in  the 
backscattered  waves  appear  to  be  bright  points.  These  bright  points  can  be  clearly  seen  in 
Figure  2.  In  a  general  radar  image,  the  coherent  summing  of  random  scatterers 
contributes  to  bright  points  on  the  imagery,  which  is  known  as  speckle  noise  [3]. 


Figure  2.  Radar  imagery  (Top)  and  a  photo  of  the  corresponding  object  (Bottom). 

(From  [4]) 

Despite  their  differences,  the  resolution  is  usually  an  important  measure  of 
performance  in  both  optical  and  radar  imagery  as  it  determines  the  ability  to  discriminate 
distinct  objects  within  the  scene.  The  factors  that  affect  radar  resolution  are  discussed 
next. 

1.  Radar  Resolution 

Radar  resolution  provides  a  measure  of  how  close  together  two  distinct  point 

targets  can  be  in  order  to  discriminate  them  as  separate.  Imaging  resolution  for  radars  is 

typically  defined  in  terms  of  the  down-range  (i.e.,  along  the  radial  direction  from  radar) 
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and  cross-range  dimensions.  This  is  illustrated  in  Figure  3.  For  a  monostatic  radar,  the 
down-range  resolution  is  dependent  on  the  pulse  bandwidth  as  ([5],  p.  13) 

ADR  =  —  (1) 

2/3 

where  / 3  is  bandwidth  of  the  transmitted  pulse  and  c  is  the  speed  of  light.  In  radars,  where 
detection  performances  are  important,  longer  pulse  widths  are  desired  as  these  improve 
detection  performances  (assuming  the  same  peak  power).  Equation  (1)  implies  that  the 
resolution  of  a  long  pulse  can  be  made  arbitrarily  small,  provided  the  pulse  is  modulated 
with  a  sufficiently  large  bandwidth  (3.  The  Linear  Frequency  Modulated  (LFM) 
waveform  is  a  common  example  of  such  modulated  waveforms  used  in  imaging. 

The  cross-range  resolution,  ACR,  for  a  stationary  monostatic  radar  depends  on  the 
angular  beamwidth,#,  of  the  radar,  which  in  turn  depends  on  the  physical  dimensions  of 
the  radar  antenna  effective  aperture  as  follows  ([6],  p.  155): 

0MB  (in  radians )  =  ~~  (2) 

where  Ka  equals  to  1 .25  for  uniform  illumination  pattern  of  a  square  antenna  aperture,  A 
is  the  carrier  wavelength  and  D  is  the  antenna’s  physical  dimension  in  the  azimuth  plane. 


Synthetic  Aperture  Radar  (SAR)  processing  allows  the  improvement  of  cross¬ 
range  resolution  of  monostatic  radars  beyond  the  bounds  of  its  physical  antenna  aperture. 
This  is  done  by  creating  a  large  synthetic  aperture  by  taking  multiple  snapshots  from  a 
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moving  platform.  This  approach  is  much  akin  to  creating  a  large  physical  planar  array, 
except  that  the  elements  are  transmitted  at  different  times.  This  scheme  also  implies  that, 
in  order  to  achieve  similar  “antenna  element  spacing,”  a  sufficiently  high  sampling  rate 
(along  the  synthetic  aperture)  should  be  needed  to  avoid  the  presence  of  artifacts  such  as 
grating  lobes  ([8],  pp.  28-29).  The  signals  obtained  along  this  “array”  are  coherently 
summed  to  achieve  the  equivalent  narrow  beamwidth.  Thus,  the  achievable  cross-range 
resolution  depends  on  the  length  of  the  synthetic  aperture  that  can  be  created.  This 
resolution  therefore  varies  with  the  different  types  of  imaging  data  collection  schemes. 
Figure  4  and  Figure  5  illustrate  two  commonly  described  schemes. 

A  strip-map  or  scan  mode  of  SAR  imaging  continually  sweeps  along  a  track  and 
produces  a  strip  of  imagery.  Here,  the  size  of  the  synthesis  aperture,  with  respect  to  a 
single  point  target,  is  limited  by  the  antenna’s  actual  beamwidth  as  it  sweeps  across  the 
track.  Hence,  the  resolution  for  a  broadside  target  is  given  by  ([5],  p.  22) 

A  CR  —  ^  ^ 

Stnpmap  2A6  2  (A/D)  2 

where  X  is  the  carrier  wavelength  and  D  is  the  physical  dimension  of  the  antenna.  This 
seems  to  be  a  counter-intuitive  result  as  it  suggests  that  the  smaller  the  antenna  dimension, 
the  better  the  cross-range  resolution.  This  is  because  a  smaller  antenna  tends  to  have  a 
wider  beamwidth,  which  results  in  a  longer  observation  time  for  a  given  target. 


Figure  4.  Stripmap  modes  for  SAR  imaging  (from  [9]) 
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In  the  spotlight  mode  of  SAR  imaging,  the  antenna  beam  is  placed  on  a  fixed 
scene  of  interest  and  the  beam  is  continually  steered  to  remain  on  this  fixed  scene  as  the 
platform  moves.  Here,  the  size  of  the  synthetic  aperture  is  typically  limited  by  the 
observation  time  and  the  radar  velocity.  Hence,  the  resolution  for  a  broadside  target  is 
given  by  [5],  p.  23. 


A  CR 
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Spotlight 


2A6  2(yT/R)  2  vT 


(4) 


where  X  is  the  carrier  wavelength,  R  is  the  range  to  the  target  in  down-range  direction,  v 
is  the  speed  of  the  SAR  platform  and  T  is  the  observation  time. 


Figure  5.  Spotlight  mode  for  SAR  imaging  (from  [9]) 


For  the  studies  done  in  this  thesis,  the  transmitters  and  receivers  are  assumed  to  be 
omnidirectional  and  range  decay  is  neglected.  Hence,  the  aperture  is  only  limited  by  the 
observation  baseline  (i.e.,  radar  velocity  times  observation  time)  simulated.  Thus,  the 
cross-range  resolution  can  also  be  given  by  Equation  (4). 

For  multi-static  radar  imaging,  the  point  spread  function  that  defines  resolution 
will  depend  on  the  multi-static  geometry  as  well  as  the  effectiveness  of  the  coherent  sum 
(e.g.,  due  to  time  synchronization).  The  down-range  and  cross-range  resolutions  are  best 
worked  out  through  simulations  [7],  [10]. 
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2.  One-Dimensional  Imaging 

In  the  simplest  case  of  a  stationary  mono-static  pulsed  radar  with  directional 
antenna,  a  single  pulse  generates  a  series  of  returns  corresponding  to  reflections  along  the 
transmit-receive  path  direction,  with  the  strength  of  the  echo  proportional  to  the  radar 
reflectivity  of  the  object  after  adjusting  for  range  decay.  This  series  of  returns  are  known 
as  a  range  profile,  and  it  provides  a  one-dimensional  image  of  the  reflectors  along  the 
down-range  axis  of  the  radar.  In  order  for  such  an  image  to  be  useful,  the  image  must 
allow  the  distinctly  separately  objects  of  interest  to  be  resolved  clearly.  Figure  6  shows  an 
example  of  a  range  profile  generated  by  a  High-Resolution  Radar  (HRR). 
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Figure  6.  Example  of  a  range  profile  generated  by  a  high  resolution  radar 

(From  [11]) 

Table  1  provides  a  feel  of  the  necessary  transmitted  pulse  bandwidth  versus 
desired  down-range  resolution.  It  is  quite  clearly  seen  that  higher  resolutions  either 
require  higher  transmission  frequencies  (and  have  higher  hardware  requirements  in  order 
to  accommodate  the  larger  bandwidths  required),  or  require  the  complex  processing  of 
waveforms  in  ultra-wideband  radar  systems,  since  narrowband  assumptions  may  no 
longer  be  valid. 
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Table  1.  Required  pulse  bandwidths  versus  desired  range  resolution  for 

mono  static  radars 


Desired  Range  Resolution 

Required  Pulse  Bandwidth 

10  meters 

15  MHz 

1  meter 

150  MHz 

0.1  meter 

1.5  GHz 

Before  proceeding,  a  brief  note  can  be  made  about  the  domain  of  validity  of  one¬ 
dimensional  (1-D)  imaging.  A  common  approximation  used  in  the  SAR  image  processing 
is  the  start-stop,  or  stop-and-go,  approximation  [12],  [13],  where  the  radar  and  target  are 
assumed  to  be  momentarily  stationary  during  their  respective  times  of  interaction  with  the 
wave.  While  literature  commonly  mentions  this  approximation  in  conjunction  with  two- 
dimensional  (2-D)  imaging,  we  note  that  the  assumption  basically  simplifies  the  1-D 
imaging  processing  (or  fast-time  domain  processing  as  it  is  commonly  known  in  SAR 
literature).  This  approximation  is  generally  valid  when  the  radar  and  target  do  not  move 
significantly  during  their  interaction  times  with  the  wave.  For  typical  radar  applications, 
the  speeds  of  the  radar  and  target  are  small  relative  to  the  speed  of  electromagnetic  wave 
propagation  and  pulse  widths  are  sufficiently  short  to  ensure  this  validity.  The 
consequence  is  that  the  received  signal  is  expected  to  just  be  a  time-shifted  scaled  replica 
of  the  transmitted  signal.  Hence,  in  the  typical  radar,  a  replica  of  the  transmitted  signal  is 
used  in  the  correlation  receiver  for  matched  filtering  to  generate  the  range  profile. 

3.  Two-Dimensional  Imaging 

Two-dimensional  radar  imaging  seeks  to  create  a  2-D  spatial  image  in  terms  of 
the  radar  reflectivity  of  an  observed  scene  or  target.  There  is  much  available  literature, 
such  as  [13],  that  provide  descriptions  of  methods  or  algorithms  to  form  a  2-D  radar 
image.  These  methods  all  seek  the  same  outcome,  which  is  to  take  advantage  of  the 
underlying  information  contained  within  observations  of  the  target  from  multiple  aspects 
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angles  to  create  an  image.  This  procedure  requires  the  observation  of  a  scene  or  target 
from  multiple  aspect  angles,  which  is  typically  achieved  using  the  following  methods: 

•  The  most  commonly  discussed  method  is  Synthetic  Aperture  Radar  (SAR) 
imaging.  This  involves  a  radar  taking  snapshots  of  the  desired  scene  at 
different  positions  while  it  moves. 

•  A  second  method  known  as  the  Inverse  SAR  or  ISAR  method  observes  a 
moving  target  with  multiple  looks  as  the  target  moves  and  turns. 

•  Multi- static  imaging  involves  having  multiple  receivers  that  receive 
scattered  waves  at  different  directions  with  respect  to  the  scene  or  target. 
The  transmitters  are  not  co-located  with  the  receivers. 

The  concept  of  “backprojection”  provides  a  simple  way  of  understanding  the 
physical  method  behind  radar  image  reconstruction.  Backprojection  essentially  involves 
taking  the  received  signal  and  mapping  it  to  the  set  of  possible  physical  locations  using  a 
known  physical  relationship.  In  the  case  of  radar,  this  is  the  relationship  between  a  spatial 
location  and  the  propagation  time  required  for  the  signal  to  travel  from  the  transmitter  to 
that  location  and  then  back  to  the  receiver.  Since  the  transmission  and  reception  times  for 
a  radar  signal  are  typically  known,  the  set  of  possible  spatial  locations  from  which  the 
associated  signal  was  backscattered  can  be  determined  1  .  The  received  signal  is 
proportional  to  the  radar  reflectivity  of  target,  weighted  by  the  beam  patterns  of  the 
transmitter  and  receiver,  propagation  effects  and  range-decay.  Hence,  the  associated  radar 
reflectivity  to  be  allocated  to  each  location  is  the  signal  strength,  adjusted  by  the  inverse 
of  the  two-way  antenna  beam  pattern  and  range-decay.  In  this  thesis,  the  antennas  are  all 
assumed  to  be  omni-directional.  Propagation  effects  and  range-decay  are  neglected  as 
only  the  relative  amplitudes  within  a  small  far-away  scene  are  examined,  and  hence  the 
differences  within  each  scene  are  small.  Hence,  the  radar  reflectivity  in  this  case  is 
directly  proportional  to  the  signal  received.  Figure  7  illustrates  an  example  of  the 
backprojected  image  for  a  single  monostatic  radar  pulse  corresponding  to  multiple  point 
targets. 


1  We  assume  that  the  radar  is  operating  within  its  unambiguous  range  and  thus  all  signal  returns  are 
from  the  pulse  that  was  transmitted  immediately  before.  In  the  case  of  a  radar  operating  with  range 
ambiguity,  the  received  signal  may  contain  a  mixture  of  returns  from  multiple  pulses  and  hence  additional 
ambiguity  is  introduced  into  the  mapping  between  the  spatial  domain  and  reception  time. 
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Figure  7.  Example  of  backprojected  image  for  a  single  monostatic  radar 
pulse  for  multiple  point  targets.  The  grey  lines  are  contours  lines  of  equal 
reception  times  passing  through  the  targets’  locations.  (From  [14]) 


To  localize  the  scatterers,  an  observation  can  be  made  from  another  aspect  angle 
and  the  backprojected  images  are  overlaid.  Figure  8  shows  the  composite  image  formed. 
It  can  clearly  be  seen  that  the  intersections  of  the  backprojections  correspond  to  the  true 
locations  of  the  scatterers.  This  is  somewhat  akin  to  position  triangulation.  In  this  case, 
the  backprojected  images  are  coherently  summed  to  obtain  the  final  image. 


Sensor  Flight  Path 


Figure  8.  Composite  backprojection  image  for  monostatic  radar  pulses  taken 
at  three  different  aspects  for  multiple  point  targets.  All  three  sets  of  lines  intersect 
at  the  true  target  locations.  (From  [14]) 


In  the  case  of  a  practical  radar  pulse,  the  response  to  each  point  scatterer  is  not  a 


delta-function,  but  rather,  spreads  over  period  of  time  around  the  target  location.  This  is 

10 


because  signals  which  are  bandlimited  in  the  frequency  domain  must  have  a  finite 
duration  in  the  time  domain.  The  consequence  of  this  in  2-D  imaging  is  that  the  imaging 
outcome  to  an  ideal  target  is  a  response  that  is  centered  on  the  target  location,  but  is 
spread  around  its  spatial  vicinity.  This  is  known  as  the  imaging  Point  Spread  Function 
(PSF).  Figure  9  shows  an  example  of  an  imaging  PSF.  The  PSF  is  of  interest  in  imaging 
studies  as  it  defines  the  imaging  characteristics  of  the  system,  including  the  imaging 
resolution  attainable.  Assuming  that  the  objects  are  non-interacting,  the  imaging  response 
to  multiple  point  targets  is  the  superposition  of  the  point  spread  functions,  each  shifted  by 
the  spatial  coordinates  and  weighted  by  the  radar  reflectivity  of  each  object. 
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Figure  9.  Example  of  an  imaging  point  spread  function  (From  [5]) 

C.  MOTIVATION  AND  FOCUS  OF  THESIS 

Radar  imaging  is  of  great  interest  in  both  civil  and  military  applications, 
especially  due  to  its  day/night  all-weather  capabilities.  The  objective  of  this  thesis  is  to 
study  the  imaging  degradations  that  occur  when  a  scene  to  be  imaged  by  radar  contains 
moving  targets.  This  work  began  as  a  study  of  the  linearized  theory  by  Cheney  and 
Borden  [12],  [15]  for  the  imaging  of  moving  targets.  During  the  course  of  investigation,  it 
was  realized  that  there  was  a  need  to  develop  an  independent  model  capable  of  modeling 
the  non-linear  effects  on  the  received  signal  resulting  from  target  and  sensor  motion,  in 
order  to  put  the  linearized  theory  to  test.  Hence,  the  focus  of  this  thesis  is  the 
development  and  verification  of  such  a  simulation  model,  in  order  to  observe  the  imaging 
degradations  that  occur  with  sensor  and  target  motion. 
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D.  THESIS  ORGANISATION 


This  thesis  is  organized  in  the  following  manner: 

Chapter  II  provides  a  theoretical  overview  of  the  physics  of  wave  interaction  with 
a  moving  target.  This  is  used  to  derive  the  signal  propagation  time,  as  a  function  of  the 
transmission  time,  for  sensors  and  targets  that  are  moving  linearly  with  constant  speed. 
The  effects  of  motion  are  briefly  explored.  This  section  concludes  with  a  short  literature 
survey  providing  an  overview  of  current  research  on  the  topic  of  radar  imaging  with 
respect  to  moving  targets. 

Chapter  III  provides  a  description  of  the  implemented  simulation  model.  The 
rationale  and  design  of  the  model  to  generate  the  received  signal  at  uniformly  spaced  time 
intervals  is  explained  in  detail.  The  implementation  of  the  backprojection  algorithm  is 
also  described. 

Chapter  IV  consists  of  two  main  portions.  First,  the  results  of  model  verification 
for  the  1-D  and  2-D  imaging  outputs  are  presented.  Then,  two  case  studies  on  the 
imaging  of  moving  targets  are  presented,  one  for  a  multi-static  configuration  and  another 
with  a  moving  SAR. 

Chapter  V  concludes  the  thesis  by  summarizing  the  results  and  provides 
suggestions  for  future  research  areas. 
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II.  MOTION  IN  RADAR  IMAGING 


In  this  chapter,  we  outline  a  mathematical  model  that  illustrates  the  scattering  of 
electromagnetic  waves  and  one-dimensional  (1-D)  imaging  of  scenes  containing  moving 
targets.  The  investigation  here  focuses  on  the  1-D  point  spread  function  (PSF) — i.e.,  the 
received  response  to  a  single  perfect  point  target  resulting  from  a  transmitted  signal.  A 
general  discussion  of  the  various  methods  found  in  literature  to  address  moving  targets 
rounds  out  the  chapter. 

A.  SCENARIO  AND  ASSUMPTIONS 

We  consider  a  mono-static  radar  system  (in  which  the  transmitter  and  receiver  are 
co-located)  mounted  on  a  moving  platform  with  trajectory  path,  (?) ,  which  in  this 
general  case  need  not  be  restricted  to  a  linear  profile.  We  also  have  a  moving  target 
whose  general  trajectory  path  is  r tgt(t)  .  To  study  complex,  extended  targets,  one 

approach  could  be  to  take  on  the  first-order  Born  approximation  which  basically  assumes 
targets  are  made  up  of  non-interacting  scatterers  (i.e.,  no  multiple  reflection  terms,  etc). 
When  such  an  approximation  is  valid,  linear  superimposition  can  be  used  to  determine 
the  returns  from  a  more  complex  target  scenario.  Beyond  this  assumption,  more 
complicated  models  that  consider  multiple -order  interactions  between  scatters  will  be 
needed  to  accurately  model  returns,  but  such  modeling  will  greatly  increase  complexity 
and  computation  requirements.  As  we  are  only  studying  a  point-spread  function,  we 
assume  a  single  perfect  point  reflector  as  our  target  in  the  scene  of  interest. 

B.  TRANSMISSION  AND  PROPAGATION 

Since  we  envision  that  the  1-D  image  is  to  be  used  as  an  input  to  a  2-D  image 
formation  process  in  SAR  processing,  we  consider  the  origin  as  the  center  of  the 
observed  scene,  within  which  the  target  exists.  At  time  t  =  0,  we  assume  the  target  and  a 

monostatic  sensor  to  be  at  position  F  (f)  =  m  (f)  and  Fte(f)  =  utx(t ) ,  respectively. 

We  assume  that  the  radar  transmits,  at  time  t  =  ttx,  an  arbitrary  single  narrowband 
radar  waveform,  w,(t),  with  pulse  width  Tw.  A  non-dispersive  transmission  medium  is 
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taken  here  as  it  simplifies  the  analysis  greatly  by  allowing  us  to  assume  that  the 
waveform  shape  does  not  naturally  change  with  time.  Without  loss  of  generality,  this 
assumption  is  equivalent  to  propagation  through  uniform  free  space,  hence  the  wave’s 
propagation  speed  is  equal  to  the  speed  of  light,  c.  From  here  on,  we  proceed  by 
neglecting  the  scaling  effects  of  range  decay  found  in  the  3-D  Green’s  function  for  point 
sources  in  free-space.  This  approximation  is  valid  when  we  are  only  looking  at  a  small 
scene  of  interest  that  is  far  away  from  the  radar,  such  that  the  scaling  effects  due  to  range 
remains  relatively  constant  within  the  scene  and  across  time  during  the  imaging  process 
([13],  p.  13).  Thus,  the  incident  EM  wave,  at  scalar  range  p  from  the  radar  and  at  time  t, 
is  given  as: 


eJpp  ;U  = 


0 


for  0<t-ttx-  —  <Tw 
c 

otherwise 


(5) 


For  simplicity  of  expression,  we  henceforth  drop  the  limits  in  time  and  just  note  that  the 
wave  goes  to  zero  outside  the  corresponding  time  arguments  defined  for  the  waveform. 

We  note  that,  when  assuming  a  perfect  omni-directional  antenna,  the  strength  of 
the  propagating  incident  wave  is  determined  by  its  range  from  the  radar  and  is 
independent  of  direction.  However,  a  scalar  range  can  no  longer  be  considered  when 
viewing  the  wave  from  the  chosen  origin,  since  the  radial  symmetry  of  propagation  is  no 
longer  valid  when  viewed  from  a  point  other  than  the  source.  As  the  radar  is  assumed  to 

be  at  r rx(ttx)  at  transmission  time,  we  can  rewrite  the  incident  wave  (i.e..  Equation  (5)) 
with  respect  to  a  position  from  origin,  p  ,  as: 


eJpPPJ  =  wr 


f  -  E 

P~  rr*(0 
t~ta - 


(6) 


This  incident  wave  is  valid  for  any  fixed  position  p  until  the  wave  reaches  an 
object  with  which  it  interacts. 
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C.  INTERACTION  AND  RECEPTION 

As  already  indicated,  for  this  study,  the  target  is  only  one  perfect  point  located  at 
YJt) .  This  target  is  also  assumed  to  be  idealized  in  such  a  way  that  it  supports  an 

infinite  bandwidth  in  its  response  and  hence  does  not  distort  the  frequency  response  of 
the  incident  wave.  It  will  thus  perfectly  and  instantly  reradiate  the  incident  wave  omni¬ 
directionally.  Thus,  the  reflectivity  is  simply  a  scaling  factor  at  the  location  at  the  target. 
For  simplicity,  we  neglect  the  negative  sign  (i.e.,  phase  shift  of  n  radians)  as  this  is  a 
constant  phase  shift  that  will  not  affect  our  results. 

Clearly,  the  target  will  interact  with  the  wave  when  the  positions  of  the  wave  and 
the  target  coincide  in  space.  For  the  sake  of  analysis,  we  denote  the  time  when  this 
happens  as  tlx,  where  1LX  >  ttx.  We  will  proceed  to  derive  the  relationship  between  this  time 
and  the  transmission  and  reception  times  for  some  specific  cases  later  in  this  chapter. 
Following  wave-target  interaction,  the  scattered  wave  that  is  re-radiated  omni¬ 
directionally  from  YT(th)  can  be  observed  at  any  given  position  q  as: 


EJq,t  ;U  =  wR 


t-t 


\q- 


-Ylx(tix) 


(7) 


The  total  field  present  in  the  environment  is  a  combination  of  the  incident  and 
scattered  waves.  However,  for  the  purpose  of  reception,  we  consider  only  the  scattered 
field.  This  is  valid  because  practical  systems  always  have  some  means  of  isolating  the 
received  signal  from  the  transmitted.  In  pulsed  systems,  the  transmit  and  listening  time 
periods  are  designed  to  be  mutually  exclusive,  and  for  CW  systems,  an  isolator  is 
normally  used.  Thus,  the  received  wave  is  just  the  scattered  wave  which  propagates  back 
to  the  receiver,  interacting  with  the  receiver  when  their  positions  coincide. 


EJYR(trx),t\t[x)  =  wR 


t  ~  - - - - - 


V 


(8) 


J 
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D.  LINEAR  TRAJECTORY 


We  note  that  in  general,  the  positions  I ~'tx(ttx)  ,  T rx(tjx)  and  (trx )  are  time 

dependent.  Hence,  the  relationship  given  in  Equation  (8)  between  the  received  wave  and 
transmitted  wave  is  a  very  general  one  and  does  not  particularly  offer  insights  into  the 
behavior  affecting  the  waveform.  However,  we  can  deduce  from  the  previous  discussion 
that  the  various  timing  within  the  argument  depends  on  the  interception  points  between 
the  wave  and  the  target/radar  and  thus  is  geometry  dependent. 

To  put  things  into  clearer  perspective,  we  select  a  case  where  both  a  monostatic 
radar  and  the  target  are  traveling  on  linear  paths  with  constant  velocities: 


r*(0  =  ua+vj 

rw(0  =  *v+v 


(9) 


Figure  10  illustrates  the  geometry  from  transmission  to  target  interaction.  As 
defined,  t-lx  is  the  time  when  the  positions  of  the  wave  and  the  target  are  equal. 


Figure  10.  Illustration  of  geometry  at  interaction  between  wave  and  target. 
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Since  we  know  the  wave  speed,  as  well  as  the  position  and  velocity  of  the  target,  the  time 
of  interaction,  tix,  can  be  computed  by  solving  for  interception,  that  is: 

R-EMwave  (  hx  ~  hx  )  —  ^Igt  (he)  (10) 


We  define  A tix  -  tlx  —  tlx and  compute  the  ranges  from  ranges  from  the  radar’s  position  at 
transmission: 


cAt-. 


U,g‘(t^)  +  ^tgAtIX-Utx(t,x) 


(11) 


For  simplicity  of  writing,  we  define  the  relative  positions  between  the  radar  and 
the  target  at  the  time  of  transmission,  i.e.,  the  expression  (utgl  +  vtgtttx  -  utx  -  vtxttx) ,  as 

it,(ftt).  Expanding  Equation  (11)  yields  the  following  equation,  which  is  quadratic  in 
terms  of  A tix : 

\U,  (fJt  +  2  ut  (; ttx )  ■  vtgtA  tix  +  (|  vtgt\-c2  )A  t*  =0  (12) 


We  solve  for  the  causal  (i.e.,  A tjx  >  0 )  solution  (See  Appendix  A. I  for  full  derivation)  and 
thus  obtain: 


UAtx)-Vtgt  +  ^ 

kut(ttx)-v,gt)2  +(c-  - 

2)|«,M2 

(c2- 

Kt 

2 

) 

(13) 


To  determine  trx,  we  apply  similar  reasoning  with  the  roles  of  the  source  and 
target  now  reversed  (i.e.,  the  target  is  now  the  source  of  the  wave  and  the  radar  is  the 

interception  point).  For  this,  we  define  A trx=trx-  tix  SO  that  trx=ttx+Atix+Atn-  We 
obtain 


At  (t,  ) 

rx  x  tx  ' 


_  (t,x  +  A 4 )  •  E,  +  >/(«,  (trx  +  A tix )  ■  vtx  f  +  (c2  -  |vft  |2 )  I u,  (t,x  +  A tix ) 


ic2-v,x  ) 


(14) 
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where  ut  (tIX  +  Atix)  can  be  obtained  by  evaluating  the  position  of  the  target  at  time  llx  =  ttx 
+  AtiX  ,  and  Atix  can  be  computed  using  Equation  (13). 

E.  WAVEFORM  DISTORTION 

From  Equations  (13)  and  (14),  we  see  that  both  linear  trajectories  (A tix  and  A tlx) 

are,  in  general,  non-linear  functions  of  ttx.  This  has  consequences  on  the  shape  of  the 
waveform.  To  visualize  this,  we  consider  a  pure  sine  waveform,  wR(t),  as  made  up  of  the 
sum  of  many  samples  separated  by  infinitesimally  small,  equally  spaced,  transmission 
time  intervals,  ttxn  =  nAttx  i.e., 


wR(t)=  lim  V  e‘27rfnMu S(t -nAt)  (15) 

Ar„  ->0  1 

n=— oo 

We  note  from  the  previous  definitions  that  the  reception  times  for  two  adjacent  samples 
transmitted  at  times  ttx  and  t,x+  Attx  can  be  expressed  as 

tM=ta+K(tJ+K(0  (i6) 

<„  + K ) = t„ + K  +  K  d„ + K  )+K«*+K)  d7) 

In  order  for  the  received  waveform  shape  to  remain  distortion-free  (i.e.,  the  same 
as  transmitted),  it  is  required  that  upon  reception,  all  samples  have  been  scaled  by  the 
same  scaling  factor  and  that  all  samples  continue  to  be  separated  by  the  time  interval,  Attx. 
This  requirement  must  be  true  for  our  free-space  propagation  and  ideal  target 
assumptions.  Thus,  for  absence  of  waveform  distortion  at  receive,  it  is  necessary  that 

Q=K  (is) 

Or,  from  Equation  (16)  and  (17),  equivalently 

K  ft,  +  K ) + K  ((« +  K) = K  ) + K  (i9) 
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We  can  see  from  Equations  (13)  and  (14)  that  these  quantities  are  non-linear 
functions  of  the  transmission  time  and  hence,  in  general.  Equation  (19)  will  not  be 
satisfied  and  we  can  expect  the  waveform  shape  to  be  distorted  due  to  changes  in  the 
spacing  of  the  samples.  This  effect  is  akin  to  time  dilation  or  expansion  due  to  relative 
motion.  Once  the  samples  are  no  longer  equally  spaced  as  per  the  original  signal,  we  can 
expect  that  wn(t)  will  no  longer  be  a  pure  sinusoidal  tone  when  viewed  in  the  frequency 

domain.  A  special  case  exists  when  both  radar  and  targets  are  stationary  ( vtx  =  0  and 

vtgt  =0),  resulting  in: 


^,-,  +  A  t, 


c 


(20) 


which  is  the  well-known  time-of-flight  formula  used  in  radar  analysis  and  is  independent 
of  the  transmission  time. 


F.  DOPPLER  MODULATION  EFFECT 

Another  common  situation  present  in  radar  calculations  is  the  scenario  where  the 
velocity  vectors  of  the  radar  and  target  are  parallel  to  their  initial  line-of-sight, 

iy0-Ftt(0  •  In  such  cases,  the  velocities  may  be  expressed  as  scalar  multiples  of 

the  initial  line-of-sight  vector.  For  simplicity,  we  assume  that  the  radar  is  stationary  and 
note  that  all  other  such  cases  can  be  expressed  as  this  basic  case  by  converting  to  the 
radar’s  frame  of  reference.  We  consider  the  case  of  an  inbound  target  so  that 

Vtt»=-H“(k)  &  (21) 


Using  Equations  (21)  in  (13)  and  (14),  we  can  see  that  the  terms  in  the  radical 
cancel  out  to  yield  a  simple  expression  (See  detailed  derivation  in  Appendix  A. Ill)  : 


K«a) 


(C  +  Kr) 


(22) 


19 


Similarly, 


At  (t,  )  = 

rx  V tx  / 


lMO 


(c+K\) 

We  can  express  the  change  in  propagation  times  between  any  two  samples  as 

Atp(ttx+nAtJ-Atp(ttx)  = 

Putting  (24)  back  into  (15),  we  have 


-2 

nK 

(c  + 

) 

Er(t)=  lim  V  cos(2nfnAt ,x)S(t-nAt  (1- 

\t.  >0  J 


2  v, 


tg>\ 


(c  +  \t) 


-)) 


=  cos(2;r/(l  + 


2  v, 


tg<\ 


(c  +  \t) 


-)t) 


(23) 


(24) 


(25) 


The  time  dilation  introduced  can  be  seen  equivalently  as  a  Doppler  frequency 
shift  of  a  radar  waveform  by  an  inbound  or  outbound  target  when  retaining  the  same  time 
spacing  at  the  receiver. 

G.  ROLE  OF  START-STOP  APPROXIMATION 


To  examine  the  effects  of  motion,  we  take  a  closer  look  at  the  implications  of  the 
commonly  assumed  start-stop  approximation  vis-a-vis  the  received  wave.  The 
fundamental  assumption  here  is  that  the  radar  and  the  target  are  both  momentarily 
stationary  when  the  wave  is  being  transmitted  and  received. 


From  the  previous  discussion,  it  is  clear  that  when  the  radar  and  the  target  are 
stationary,  the  received  signal  will  be 


E,-x (f  r (trx ),t)  =  wR(t- tlx  - tp )  for  tlx  <  t  <  tlx  +  Tw 


where  tp  -  2 


rr(0-r«(0 


/  c 


(26) 


This  means  that  at  the  receiver,  the  correlator  can  be  implemented  by  performing 
matched  filtering  with  a  replica  of  the  transmitted  waveform.  The  output  of  the  matched 
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filter  is  given  by  the  autocorrelation  function  of  WR(t) ,  given  by  Equation  (27),  which 


has  its  peak  when  r  =  ttx  + 1  . 

p(j)  =  \wR(t-tlx- tp )  w*  (f  - t)  dt  (27) 

where  the  asterisk  denotes  complex  conjugation. 

When  the  radar  and  target  are  moving,  the  distorted  waveform  will  change  the 
output  of  the  correlation  receiver.  The  exact  effect  of  the  distortion  depends  on  the  type 
of  waveform  employed.  This  is  typically  studied  through  the  use  of  the  radar  ambiguity 
function,  which  is  essentially  the  correlation  function  in  the  spatial  domains,  evaluated  at 
different  modulation  frequencies,  v,  as  follows 

oo 

X(v,  t)=\wr  ( t)wR  ( t  -  t)  exp(2 nvt)dt  (28) 


H.  ARTIFACTS  FROM  MOVING  TARGETS 

Imaging  degradations,  also  termed  artifacts,  occur  when  the  conditions  of  the 
actual  scene  do  not  correspond  to  the  assumptions  embodied  within  processing 
techniques.  One  such  mismatch  is  the  assumption  of  a  stationary  scene  when  there  is  in 
fact  motion  present.  Typically,  the  moving  radar  will  have  its  own  means  of  knowing  and 
compensating  for  its  own  motion,  for  example  through  the  use  of  Inertial  Navigation 
Systems  (INS)  and/or  Global  Positioning  System  (GPS),  or  through  antenna  phase  center 
displacement.  Hence,  available  literature  mainly  focuses  on  addressing  the  target  motion 
and,  to  a  lesser  extent,  the  residual  effects  of  imperfect  radar  motion  compensation.  We 
focus  solely  on  target  motion  here. 

Qualitatively,  the  outcome  of  target  motion  in  a  2-D  SAR  image  is  generally  as 
shown  in  Figure  11.  These  generally  include  smearing  in  range,  smearing  in  cross-range, 
and  a  shift  in  the  cross-range  dimension. 

The  analysis  from  [16]  and  ([17],  pp. 152-154)  illustrates  that  the  effects  of  the 
waveform  distortion  leads  to  the  shift  in  the  cross-range  dimension  due  to  the  Doppler 
shift  introduced  and  a  cross-range  defocusing  due  to  a  non-linear  phase  term.  Other 
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additional  effects  include  the  effect  of  range-walk  (i.e.,  target  crossing  into  a  different 
range  cell  during  the  processing  interval)  and  the  effect  of  observing  the  target  at 
different  positions  over  the  observation  interval  in  the  backprojection  algorithm.  We 
attempt  to  study  these  effects  using  the  simulation  model  developed. 


SAR  image  of  moving  target 


Defocused  SAR  image 
of  a  moving  target 


Target's  true 
location 

© 

T 


Smeared 
in  range 


Smeared 

in  cross  range 
t 
i 


Shifted 


0» 

3 

CD 

n> 


Cross  range 

Figure  11.  Effects  of  target  motion  on  SAR  imaging  (From  [17]) 


I.  METHODS  TO  ADDRESS  DE-FOCUSING 

The  topic  of  moving  targets  has  been  the  focus  of  ongoing  research  interest  on 
SAR  over  the  past  decades.  A  general  survey  of  literature  was  made  to  provide  brief 
backgrounds  on  past  and  current  work  in  this  area.  Based  on  this  survey,  related  research 
can  be  generally  classified  into  the  following  areas. 

1.  Velocity  Estimation  Techniques  in  Single  Channel  SARs 

Such  SAR  systems  have  only  a  single  transmit/receive  channel,  and  hence  the 
Doppler  shifted  signals  from  moving  targets  are  mixed  in  with  the  returns  from  stationary 
targets  as  part  of  the  received  signal.  During  matched  filtering  (typically  matched  for 
stationary  scene),  the  returns  from  moving  targets  will  lead  to  artifacts  as  a  result  of  their 


22 


mismatch  with  the  filter.  With  such  systems,  the  focus  is  thus  on  obtaining  an  estimate  of 
the  target  velocity  using  the  radar  returns  [18],  [19],  [20]  before  applying  the  estimated 
velocity  to  achieve  Doppler  correction  for  the  moving  target.  While  the  moving  target 
becomes  focused  with  such  corrections,  the  rest  of  the  scene  not  at  the  same  velocity 
tends  to  become  defocused. 

2.  Multichannel  SAR  (MSAR) 

This  class  of  SARs  can  generally  be  considered  to  include  the  Velocity  SAR  or 
VSAR  proposed  in  [21],  Displaced  Phase  Center  Array  (DPCA),  as  well  as  further  work 
on  the  Multichannel  SAR  (MSAR)  utilising  Space  Time-Adaptive  Processing  techniques 
(STAP)  in  SAR  [22].  The  key  characteristic  of  this  class  of  techniques  is  that  the  SAR 
system  is  set  up  to  receive  signals  in  multiple  elements  of  the  antenna  arrays  (space 
samples)  and  over  different  instances  (time  samples).  Farina  and  Lombardo  [22]  provide 
a  good  overview  of  the  evolution  of  MSAR  architectures  all  the  way  through  to  the 
eventual  incorporation  of  STAP  techniques,  including  the  use  of  time-frequency  analysis 
methods  [17].  MSARs  are  reported  to  provide  better  clutter  suppression  due  to  the 
additional  degrees  of  freedom  offered  by  the  multiple  antennas.  This  lends  itself  to 
Moving  Target  Indicator  (MTI)  processing;  however,  this  is  gained  at  the  expense  of 
greater  hardware  requirements. 

3.  Compressive  Sensing 

Compressive  Sensing  (CS)  seems  to  be  a  recently  explored  technique  introduced 
to  the  domain  of  SAR  processing  [23],  [24],  [25].  According  to  Patel  [26],  “compressed 
sensing  enables  the  reconstruction  of  sparse  or  compressible  signals  from  a  small  set  of 
non-adaptive,  linear  measurements.  If  properly  chosen,  the  number  of  measurements  can 
be  much  smaller  than  the  number  of  Nyquist  rate  samples.”  Essentially,  CS  involves  the 
selection  of  an  appropriate  sub-space  basis  set  of  functions  that  allow  the  SAR  data  set  to 
be  reduced.  This  method  has  been  applied  to  the  SAR  imaging  of  uniformly  moving 
targets  in  [23]  “with  higher  resolution  and  lower  sidelobes  with  fewer  measurements.” 
We  will  not  further  investigate  the  application  of  CS  in  thesis. 
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4. 


Linearized  Theory  by  Cheney  and  Borden 

The  linearized  theory  from  Cheney  and  Borden  [15]  is  of  particular  interest  to  this 
thesis,  as  this  thesis  was  motivated  by  study  of  this  theory.  This  imaging  procedure  for 
moving  targets  is  derived  from  the  principles  of  wave  scattering  that  combines  the  spatial, 
temporal  and  spectral  aspects  of  scattered  field  based  imaging.  The  work  is  an  extension 
of  the  prior  work  from  [12],  which  considered  stationary  sensors  in  multistatic 
configurations,  to  moving  sensors.  In  essence,  a  linearized  formulation  of  the  phase  (i.e., 
time  argument)  of  the  received  signal  is  developed  by  making  the  following  assumptions: 

•  The  Born  approximation  is  used  to  express  the  scattered  wave  as  a  simple 
function  of  the  incident  wave,  since  scatters  are  assumed  to  be  non¬ 
interacting. 

•  The  Platform  and  target  motion  are  slow  relative  to  the  speed  of  light.  The 
sensor  is  also  much  more  distant  than  target  from  the  origin  and  this  range 
is  much  larger  that  the  distance  travelled  between  pulses  by  the  sensor  or 
the  target.  This  approximation  allows  the  higher  order  terms  in  the  Taylor 
series  expansion  of  the  time  dependent  range  between  target  and  sensor  to 
be  dropped,  thus  linearizing  the  argument  that  appears  in  the  phase 
function  of  the  received  wave. 

With  the  following  definitions  adopted  in  [15], 

fm(t  -T*)  is  the  waveform  for  m-th  pulse  sent  out  at  time  t  =  Tjn 

ym ,  ym  :  transmitter  position  and  velocity  vectors  at  the  time  of  transmission 
Ym ,  :  transmitter  location  and  velocity  vectors  at  the  time  of  reception 

rm  +  vj  is  the  trajectory  vector  of  the  target 

Q( z)  is  the  3-D  spatial  reflectivity  function  of  target 

Tm  is  the  time  of  where  the  transmitted  wave  interacts  with  the  target 
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These  assumptions  lead  to  the  formulation  of  the  imaging  algorithm  for  each  pulse  based 
on  matched  filtering  that  is  based  on  the  signal  that  is  expected  from  each  point  p  and 

each  velocity  vm  in  which  the  image  Im  is  given  by 


SP,vm)  =  \OASt,  pS’J)(A7i)2 


K 


Yn. 


+  (29) 


where 


1  +  0, 


T  k,„  ([i  -  p:  v  -  k  (z,  vm  >  /  c>  -  vm )  /  o]  -  t 


riT  •  T  /  n  nR  *R  / 

Pm  =rm-rm/c  &  Pm  =rm-ym/c 


(l  — 

v„,,m  t  ,  *r 


1  +  yT  ■  v  /  c 

i  m  m 


1  +  7  -V  /c 

/  m  m 


rr(v  v  )=  vR -9r-(z+t  -vT'+yrtr) 

m  m  /  /  m  /  m  ^  m  mm  /mm' 


R1  (z,v  )=  /  .(z  +  r  -v  T  +  v  T  ) 

m  V  7  m 7  /  m  /  m  v  m  mm  *  m  m  / 


(30) 

(31) 

(32) 

(33) 

(34) 


This  method  is  similar  to  the  “shift-varying  filter”  mentioned  in  ([13],  pp.  212-215),  but 
further  accounts  for  shifts  in  the  velocity  domain,  rather  than  just  spatial  locations.  Full 
treatment  of  the  derivation  for  the  Cheney-Borden  procedure  can  be  found  in  [15].  The 
earlier  work  in  [12]  serves  as  a  useful  reference  to  provide  further  insights  to  the  basis  for 
the  derivation. 
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III.  DESCRIPTION  OF  MODEL 


This  chapter  provides  design  considerations  and  descriptions  of  the  model  that 
was  implemented  for  the  purpose  of  these  studies.  The  model  was  implemented  using 
MATLAB,  R201  la  32-bit  version. 

A.  OVERVIEW  OF  SINGLE  PULSE  MODEL 

A  key  objective  of  this  model  is  to  allow  the  synthesis  of  a  received  waveform 
that  has  been  distorted  by  effects  of  motion  for  analysis  of  effects  on  radar  images. 
Figure  12  illustrates  the  general  concept  behind  the  model  for  single  pulse  generation. 


Range  Profile 


Figure  12.  Block  diagram  for  range  profile  generation 
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B.  DEFINITION  OF  TARGET  AND  SENSOR  TRAJECTORY 


The  implemented  model  allows  for  the  definition  of  position  and  motion  in  three- 
dimension  (3-D)  spatial  coordinates.  As  this  model  is  intended  for  imaging,  the  origin  is 
set  at  the  center  of  the  imaged  scene.  This  model  will  only  work  for  static  or  constant 
linear  velocity  scenarios  due  to  the  assumption  made  during  the  derivation  of  wave- 
interaction  time  and  reception  time.  However,  the  positions  of  the  transmitter,  target,  and 
receiver  can  be  defined  independently.  This  design  is  intended  to  allow  for  the  possibility 
of  defining  different  sensor  configurations  in  different  studies.  Possible  implementation 
methods  for  different  configurations  are  provided  in  Table  2. 


Table  2.  Examples  of  possible  sensor  configuration  that  can  be  studied 


Configuration 

Possible  Implementation  Method 

Stationary,  Mono-static 

Transmitter  and  receiver  initial  positions  are  defined  to  be 

the  same,  both  with  zero  velocity.  Each  run  provides  the 

range  profile  for  one  pulse.  This  configuration  may  not  be 

useful2  for  SAR  studies  by  itself,  but  serves  as  a  useful 

case  for  model  verification  or  a  study  on  the  performance 

of  networked  monostatic  sensors. 

Stationary,  Multi-static 

Transmitter  and  receiver  initial  positions  are  defined  to  be 

different,  both  with  zero  velocity.  This  can  be  used  for 

multi-static  SAR.  In  case  of  multiple  receivers,  a  single 

pulse  for  each  receiver  position  can  be  simulated  separately 

for  the  same  transmission  instant  to  obtain  the  set  of  range 

profiles  for  a  single  image 

-  Although  such  a  configuration  can  be  practically  used  for  Inverse-Synthetic  Aperture  Radar  (ISAR), 
it  is  not  useful  here  since  the  constraint  of  linear  motion  limits  target  rotational  motion  that  ISAR  requires. 
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Moving,  Single-Channel 

Monostatic 

Transmitter  and  receiver  initial  position  are  defined  to  be 

the  same,  both  with  the  same  non-zero  velocity.  This  can 

be  used  for  single  channel  SAR. 

Moving,  Multi-channel 

Transmitter  and  receiver  initial  position  are  defined  to  be 

Monostatic 

separated  by  the  antenna  displacement,  both  with  the  same 

(e.g.,  Displaced  Phase 

non-zero  velocity.  Single  pulse  for  each  receiver  position 

Center  Arrays) 

to  be  simulated  separately  for  the  same  transmission  instant 

to  obtain  the  entire  set  of  range  returns  for  a  single  pulse. 

This  can  be  used  for  multi-channel  SAR. 

C.  WAVEFORM  GENERATION  CONCEPT 

In  order  for  the  general  waveform  model  to  work,  the  continuous-time  function  of 
the  waveform  is  required  to  be  implemented  so  that  the  complex  amplitude  of  the  wave  at 
any  point  in  time  can  be  computed.  The  reason  for  this  will  be  clearer  by  the  end  of  this 
section.  These  formulations  are  generally  well  known  (e.g.,  in  [27])  for  typical  practical 
waveforms  for  radar  systems,  and  hence  should  typically  not  present  difficulty  in 
implementation.  Examples  for  the  rectangular  and  Linear  Frequency  Modulation  (LFM 
or  commonly  known  as  chirp)  waveforms  are  provided  in  Table  3. 


Table  3.  Examples  of  waveform  functions 


Name  of  waveform 

Expression 

Rectangular  Pulse 

wR  (0  =  Arect( -  ^)  exp(i2nfct) 

vv 

Linear  Frequency  Modulated 

(Up-chirp) 

wR(t)  =  Arect(  )exp(  in  w  t  )exp(i2nfct) 

w  w 
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We  also  note  that  in  general,  the  transmitted  and  received  waveforms  are 
continuous  time  functions,  but  for  the  purpose  of  digital  simulations,  these  can  only 
modeled  through  discrete  time  functions,  typically  by  using  evenly-spaced  time  samples. 
The  use  of  evenly-spaced  time  samples  is  convenient  as  this  allows  easy  signal  and 
spectral  analysis  using  the  Discrete  Fourier  Transform  (DFT)  [28].  This  simplification 
results  in  more  cost-efficient  hardware  implementation.  As  a  result,  the  Digital-to-Analog 
Convertors  (DAC),  Analog-to-Digital  Convertors  (ADC),  and  signal  processing  units 
employed  as  part  of  the  sensor  hardware  largely  use  evenly-spaced  time  samples. 

We  see  from  the  theory  covered  in  Chapter  II  that,  for  moving  sensors  and  targets, 
the  propagation  time  for  a  wave  is  dependent  on  the  geometry  at  the  instants  of  wave 
transmission,  interaction,  and  reception.  For  constant  linear  sensor  and  targets  velocities, 
the  relationship  between  the  propagation  time,  Atp,  for  each  sample  can  be  determined  as 
the  sum  of  Equations  (13)  and  (14).  This  propagation  time,  Atp,  is  in  general  a  non-linear 
function  of  the  transmission  time.  The  implication  of  this  non-linearity  is  that,  following 
propagation,  the  corresponding  amplitude  samples  arriving  at  the  receiver  may  no  longer 
be  spaced  at  the  same  time  intervals  as  the  transmitted  waveform,  and  may  even  be  non- 
evenly  spaced  in  time  with  respect  to  each  other.  Figure  13  illustrates  an  example  of  the 
relationship  between  the  transmitted  waveform  samples  and  the  distorted  received 
waveform  samples. 

In  digital  correlation  reception,  the  cross-correlation  is  used  to  measure  the 
similarity  between  the  sampled  received  signal,  wrx(n),  and  the  matched  filtered ,  fm(n), 
and  is  computed  as: 

oo 

Prx,Jl)  =  X  wrx(n)fm(n~ 0  forleZ  (35) 

n=— co 

where  n,  l  indicate  sample  sequence  numbers  for  each  signal  sequence.  The  matched 
filter,  fm(n),  is  typically  a  sampled  replica  of  the  transmitted  waveform. 
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Figure  13.  Illustration  of  relationship  between  transmitted  and  received 
waveforms 

It  is  important  to  note  that  the  summation  in  Equation  (35)  makes  the  assumption 
of  equally  spaced  samples  in  both  the  received  and  matched  filter  sequences.  When  the 
received  samples  do  not  arrive  with  the  same  spacing  in  time  as  the  transmitted  waveform, 
these  samples  cannot  be  used  directly  in  the  correlation  receiver  to  derive  a  meaningful 
range  profile,  since  the  transmitted  and  received  time  samples  in  fact  exist  at  different 
instants  in  time,  and  thus  do  not  contribute  to  the  correlation  summation  correctly  if 
equation  (35)  is  simply  applied.  Hence,  it  is  clear  that  some  book-keeping  in  the  time 
domain  needs  to  be  done,  in  order  for  the  received  signal  to  be  processed  correctly  in  a 
discrete  simulation  of  the  received  waveform.  Based  on  the  preceding  discussion,  there  is 
therefore  good  motivation  to  consider  “resampling”  the  received  signals  to  evenly  spaced 
time  intervals  that  match  those  of  the  transmitted  waveform. 

A  number  of  methods  for  this  “resampling”  were  initially  considered.  One 

possibility  is  to  interpolate  the  received  amplitude-time  samples  to  obtain  the  required 

evenly  spaced  values  of  the  received  wave.  In  theory,  sine-interpolation  provides  a 

unique  and  perfect  reconstruction  of  the  bandlimited  signals  which  have  been  sampled  (at 

least)  at  Nyquist  rate.  The  time  shift  properties  of  the  Discrete  Fourier  Transform  (as 
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given  in  Equation  (36))  and  “spectral  interpolation  (i.e.,  zero-padding  in  the  frequency 
domain)  also  provide  convenient  means  to  obtain  the  value  of  a  function  at  a  different 
time  instant  through  a  given  sampled  time  sequence. 

x(/2-/),(mod/V) ,  DFT  sX(k)ei2*k'  (36) 

where  N  is  the  length  of  the  sequence,  I  is  the  number  of  samples  to  be  shifted  by  and  n,k 
refer  to  the  indices  of  the  processed  sample.  These  methods  do  not  work  when  the 
arriving  samples  are  non-uniform  in  time.  The  problem  of  interpolation  of  signals 
arriving  at  non-uniform  instances  in  time  is  still  an  active  area  of  ongoing  research  [28], 
[29] .  Consequently,  we  do  not  consider  these  methods  for  implementation. 

Another  possibility  is  to  make  use  of  curve-fitting  to  achieve  the  interpolation. 
The  values  of  ttXi„  and  trx,n  are  known.  Hence,  it  could  be  possible  way  to  consider  trxn  as 
an  arbitrary  function  of  trXjtl  and  search  for  a  function  that  fits  the  known  data  points.  The 
forms  Equations  (13)  and  (14)  suggest  that  the  functional  relationship  is  likely  to  include 
non-integer  powers  of  ttx,n ■  While  it  is  possible  to  reconstruct  the  curve  if  the  right 
parameters  can  be  found  for  this  curve  fitting,  the  search  space  for  these  parameters  is 
likely  be  large  due  to  the  degree  of  freedom.  Hence,  this  approach,  although  theoretically 
possible,  may  be  impractical. 

Instead,  we  choose  an  approach  based  on  the  physical  nature  of  the  problem.  For 
constant  linear  velocities,  we  can  see  that  the  geometry  between  transmitter,  the  target 
and  the  receiver  is  defined  at  the  instant  of  transmission.  By  including  the  constraint  of 
causality,  we  arrive  at  a  unique  solution  mapping  the  transmission  time,  ttXill,  to  the 
reception  time,  trXln  through  the  mathematics  laid  out  in  Chapter  II.  We  now  assert  that  a 
wave  transmitted  at  a  later  instant  will  not  be  received  earlier  than  the  earlier  transmission 
for  a  fixed  geometrical  relationship.  This  can  be  easily  shown  by  noting  that,  for  any 
given  target  speed  v,  the  propagation  time  will  decrease  most  rapidly  from  one 
transmission  instant  to  the  next  when  the  target  is  moving  directly  towards  the 
transmission  source.  In  this  scenario,  assuming  the  initial  range  separation  between 
source  and  target  is  R ,  we  note  from  Equation  (22)  and  (23)  that  the  two-way  propagation 
time  is 
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2  R 


(37) 


c  +  v 

Hence,  the  reception  times  due  to  transmission  at  two  instances  separated  by  a  small  time 
interval  Attx  are 


Kx\  hx ' 


2  R 


c  +  v 


Kx  2  -  ftx  +  + 


R  -  vA  t 


c  +  v 


Kxl  Kxl  ^tx  ' 


vAt.. 


■  =  Ac 


^  c  ^ 


c  +  v 


>0 


vc  +  v; 


V  At(A  >  0 


(38) 

(39) 

(40) 


The  difference  in  reception  times  in  Equation  (40)  is  always  greater  than  zero  for 
positive  values  of  A tn ,  hence  proving  that  reception  time  is  a  monotonic  increasing 

function  of  transmission  time.  With  these  observations,  we  can  formulate  a  simple 
iterative  search  algorithm  to  find  the  values  of  transmission  times  that  will  provide 
evenly-spaced  reception  times  as  follows: 


For  each  value  of  reception  time  to  computed 

Step  1:  Set  the  desired  reception  time  to  be  Trx  n  =  trx  ()  +  nAtlx  where  trXto  is  the 

reception  time  of  the  first  sample,  n  is  an  integer  and  Attx  is  the  sampling  interval 
used  to  generate  the  transmitted  wave. 

Step  2:  Compute,  as  an  initial  reception  time,fn  n,  using  equations  (13)  and  (14) 
with  the  initial  estimated  transmission  time  tCurrer"  =  ttx  0  +  nAttx . 

Step  3:  Compare  the  estimated  reception  time,  tnn,  against  the  desired  reception 
time,  and  do  the  following: 
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•  If  yrx,n  ~  Kx,n  \  >  \  ,  where  Ae  is  the  acceptable  error,  set  the  updated 
estimated  transmission  time  as 

f New  Current  _J  x 

txn  tx  n  v  rx,n  rx,n ' 

Compute  the  updated  reception  time,? ,  using  equations  (13)  and 
(14)  with  the  updated  estimated  transmission  time  tNew  and  repeat 
Step  3. 

•  Else,  end  the  iteration.  tCurrent  is  the  required  value  for  the  current 
sample. 

Step  4:  Repeat  from  Step  1  for  next  sample  until  all  samples  in  the  range  profile 

have  been  processed. 

Once  the  entire  set  of  corresponding  transmission  times  have  been  determined,  the 
functional  relationship  of  the  transmitted  waveform  can  be  used  to  obtain  the  associated 
amplitude  values.  It  is  for  this  reason  that  the  continuous  time  representation  of  the 
waveform  is  a  pre-requisite  for  this  method  to  work.  Figure  14  illustrates  the 
implemented  process. 

We  also  note  that  in  conventional  implementations  of  digital  matched  filtering, 
the  down-converted  received  signal  is  generally  sampled  at  uniform  rates.  Hence,  the 
received  waveform  generated  would  also  be  an  approximate  representation  of  waveforms 
received  in  actual  systems;  however,  in  practice,  there  are  other  practical  effects  that  will 
affect  the  received  signal. 
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Figure  14.  Illustration  of  received  waveform  generation  process 
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D.  MATCHED  FILTERBANK  AND  RANGE  PROFILE 


The  intent  of  matched  filtering  is  to  maximize  the  Signal-to-Noise  Ratio  (SNR) 
ratio  at  the  receiver  end  in  the  presence  of  noise.  When  only  white  noise  is  considered,  it 
can  be  shown  ([30],  pp.  22-28)  that  the  ideal  matched  filtered  is  a  template  of  received 
signal  free  from  white  noise,  but  including  distortions  such  as  Doppler  effect  encountered 
during  propagation.  Matched  filtering  in  radar  is  commonly  understood  as  using  a  replica 
of  the  transmitted  signal  or  its  Doppler-shifted  version  to  recover  the  transmitted  signal 
amidst  noise.  However,  this  model  is  not  always  true,  although  it  is  frequently  assumed  in 
radar  signal  processing.  This  matched  filter  is  a  replica  of  the  transmitted  signal  or  its 
Doppler-shifted  version  only  in  scenarios  with  stationary,  radially  ingressing,  or 
egressing  targets  when  viewed  from  the  sensor’s  frame  of  reference.  More  generally,  the 
formulation  in  [15]  provides  a  linearized  model  that  can  be  used  to  derive  the  matched 
filter  that  accounts  for  phase- space  (or  equivalently,  time  domain)  distortions  due  to 
target  and  sensor  motions.  Nevertheless,  many  systems  continue  to  use  the  transmitted 
signal  replica  or  its  Doppler-shifted  version  as  this  greatly  simplifies  signal  processing 
implementations,  yet  at  the  same  time  yields  relatively  accurate  results. 

In  this  study,  as  the  intent  is  to  examine  the  artifacts  due  to  motion,  we  only 
implement  the  matched  filter  as  a  bank  of  the  transmitted  signal  replica  and  its  Doppler- 
shifted  versions  to  observe  the  outcomes  of  wrongly  assuming  the  start-stop 
approximation.  With  such  a  bank  of  Doppler- shifted  filters,  the  output  of  the  correlation 
filter  in  the  absence  of  distortion  is  simply  the  waveform’s  ambiguity  function.  This  also 
provides  a  useful  check  for  the  correctness  of  the  simulation.  To  provide  the  best  signal 
for  imaging,  the  output  from  the  filter  within  the  banks  that  provides  highest  peak  should 
be  chosen  as  the  range  profile  to  be  used  for  imaging  since  this  filter  is  the  one  that  best 
corrects  for  the  signal  modulation  that  has  been  introduced  into  the  received  signal.  We 
note  that  the  simulation  is  sufficiently  flexible  to  allow  the  use  of  other  forms  of  matched 
or  mismatched  filters  that  the  user  generates. 

The  output  of  the  digital  correlation  receiver  can  be  described  by  Equation  (35). 
In  a  practical  simulation,  suitable  limits  need  to  be  placed  on  the  time  domain  for 

summation.  For  calculating  the  cross-correlation  between  two  sequences  of  lengths  M 
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and  N  respectively,  the  final  sequence  length  will  be  M+N-l.  For  a  time-limited 
waveform,  we  can  make  use  of  the  fact  that  the  signal  is  zero  outside  the  pulse  width. 
Hence,  the  matched  filter  length  is  determined  by  pulse  width  divided  by  the  sampling 
interval.  The  length  of  the  received  waveform  will  need  to  consider  the  time  expansion  or 
compression  introduced  during  propagation.  We  estimated  the  required  length  by  taking 
the  worst  case  time  expansion  to  ensure  that  the  entire  received  signal  is  modeled.  This 
can  be  done  by  deriving  the  time  expansion  for  the  case  of  an  egressing  target.  We  obtain 

vA  t  (  c  ) 

Kxi-Kxi  =At,x  + - ~  =  Attx  -  >0  V  Ata  >  0  for  v  <c  (41) 

c-v  \c-vj 

Hence,  for  a  transmitted  signal  sequence  with  M  samples  within  its  pulse  width,  at  least 
M  (c  /  c  -  v)  samples  (including  zero-valued  samples  outside  the  pulse  width)  should  be 
generated  for  the  received  pulse,  in  order  for  the  received  waveform  to  be  completely 
simulated.  We  add  to  this  discussion  a  few  words  about  the  importance  of  choosing  the 
right  sampling  rate.  In  general,  in  order  to  reproduce  a  signal  without  aliasing,  the 
required  sampling  rate  should  be  at  least  twice  the  signal  bandwidth  (i.e.,  Nyquist  rate). 
This  choice  of  sampling  needs  to  consider  the  added  bandwidth  due  to  the  effects  of 
distortion  in  order  to  avoid  aliasing.  In  our  study,  we  select  a  sampling  rate  much  higher 
than  the  Nyquist  rate  to  avoid  this. 

E.  IMAGE  RECONSTRUCTION 

In  image  reconstruction,  the  objective  is  to  create  a  representation  of  the 
reflectivity  of  the  scene  in  the  2-D  spatial  domain  by  using  several  observations  of  the 
desired  scene.  So  far,  the  generation  of  range  profile  outputs  needed  for  image 
reconstruction  have  been  covered  in  the  previous  sections.  Image  reconstruction  is 
discussed  here.  The  concept  of  image  reconstruction  is  presented  in  Figure  15  and 
requires  two  important  processes,  namely,  Range  Alignment,  and  Image  Formation. 
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Figure  15.  Block  diagram  illustrating  concept  of  image  reconstruction 


1.  Range  Alignment 

Typically  in  radar  imaging,  it  is  necessary  to  range  align  all  the  range  profiles 
prior  to  passing  the  data  set  to  the  algorithm.  This  is  required  as  the  radar  data  is  collected 
at  different  start  ranges  and  hence  the  same  sample  numbers  in  each  of  the  range  profile 
sequences  correspond  to  different  reception  times  (and  hence  spatial  samples).  The  intent 
of  range  alignment  is  to  ensure  that  when  processed,  the  respective  samples  in  the  range 
profiles  represent  the  same  range  when  processed  together  side-by-side.  In  order  to  range 
align  the  data,  the  amount  of  range  shift  required  can  be  worked  out  by  keeping  track  of 
the  propagation  time  to  the  scene  center  and  subtracting  this  from  the  reception  time  of 
the  signal  from  the  target.  When  this  is  done,  what  remains  is  the  time  position  of  each 
sample  relative  to  the  scene  center.  Thus,  the  signal  sample  with  time  reference  nearest  to 
zero  will  be  the  nearest  sample  to  the  scene  center  in  terms  of  range  (or  bistatic  range  for 
multi-static  radars).  All  the  different  range  profiles  can  thus  be  aligned  by  “placing”  this 
sample  at  the  middle  of  each  sequence  and  zero-padding  the  sequence  where  necessary  in 
order  to  achieve  the  same  sequence  length  for  all  the  range  profiles  to  be  processed.  As 
the  range  profiles  consist  of  evenly-spaced  samples,  the  range  alignment  to  correct  for 
fractional  steps  (i.e.,  to  shift  by  fractions  of  a  single  range  cell)  can  be  done  by  using  the 
Fourier  time  shift  property  presented  in  Equation  (36).  We  note  that  this  range  alignment 
is  exact  for  stationary  sensors  (mono-static  and  multi-static)  since  the  time  to  scene  center 
remains  unchanged  through  all  the  collected  samples.  This  is  no  longer  true  in  the  case 
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with  moving  sensors;  in  such  cases,  the  time  to  scene  center  at  the  start  of  transmission  is 
used  as  an  approximation.  We  note  that  a  separate  range  alignment  process  is  no  longer 
required  when  using  the  implemented  backprojection  method  for  image  formation. 

2.  Image  Reconstruction  Using  Backprojection 

The  concept  of  backprojection  was  presented  in  Chapter  I.  One  implementation  of 
a  backprojection  algorithm  is  the  Inverse  Radon  Transform  that  is  implemented  in 
MATLAB  as  the  function  iradon.m  [31].  The  Inverse  Radon  Transform  is  used  in 
computerised  tomography  and  is  equivalent  to  the  backprojection  of  a  scene  to  form  a 
spatial  image  when  given  a  series  of  one-dimensional  views.  This  is  similar  in  application 
to  the  use  of  a  series  of  range  profiles  that  are  collected  at  different  aspect  angle  to  a 
target,  in  order  to  form  an  image  of  the  target  reflectivity  distribution. 

A  further  simplification  of  the  Inverse  Radon  Transform  involves  using  the  two- 
dimensional  Inverse  Discrete  Fourier  Transform  (2-D  IDFT)  to  achieve  imaging.  This 
method  can  be  shown  to  be  approximately  equal  to  the  back  projection  method  under  the 
assumption  that  the  observation  angle  is  small.  This  is  because  when  the  observation 
angle  is  small,  the  data  collected  appears  on  a  nearly-rectangular  grid,  and  hence,  is 
suitable  for  processing  using  Fourier  methods  as  a  spatial  transform.  When  the  data 
collection  deviates  sufficiently  from  a  rectangular  grid  (monostatic  radar  data  are 
normally  collected  on  a  polar  grid  and  multi-static  data  on  an  elliptical  grid),  pre¬ 
processing  to  interpolate  the  data  to  a  rectangular  grid  should  be  done  prior  to  using  the 
FFT  method  for  imaging. 

We  note  that  the  Inverse  Radon  Transform  implemented  in  MATLAB  is  limited 
to  the  use  of  uniformly  spaced  observation  aspect  angles.  The  use  of  2D-FFT  also 
requires  evenly  spaced  sampling  in  both  the  range  and  angular  domains.  Since  it  is 
impractical  to  expect  radar  data  to  be  always  collected  from  uniformly  spaced  aspect 
angles  in  both  multistatic  and  moving  SAR  configurations,  an  imaging  algorithm  was 
developed  and  implemented  as  part  of  this  study  to  overcome  this  limitation. 

The  key  concept  behind  this  imaging  algorithm  is  also  the  backprojection  method. 
This  algorithm  is  inspired  by  existing  literature  describing  this  method  ([13],  pp.357- 
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359).  However,  the  implementation  here  was  developed  and  tailored  to  fit  the  rest  of  the 
model.  The  rationale  and  implementation  of  this  model  are  further  explained  here. 

Physically,  the  received  signal  is  the  result  of  the  transmitted  pulse  travelling  to  an 
object  and  being  backscattered  to  the  receiver.  The  reception  time  is  simply  the 
transmission  time  with  two-way  propagation  time  of  the  receiver.  For  each  signal 
received  at  a  given  reception  time,  the  propagation  time  can  be  deduced  if  there  is 
knowledge  of  the  transmission  time.  It  is  therefore  possible  to  reconstruct,  based  on  the 
transmitter  and  receiver  geometry,  the  set  of  possible  spatial  locations  that  this  signal  has 
been  backscattered  from. 

For  a  monostatic  radar,  this  procedure  involves  solving  for  the  set  of  possible 
position  coordinates  [x,y,z]  that  satisfy  the  following: 

tp  =  2<J(x-xr)2 +(y-yr)2 +(z-zr)2 I  c  (42) 

where  [xr,  y,;  zr]  are  the  position  coordinates  of  the  radar.  This  turns  out  to  be  the 
equation  for  circles  with  the  radar  at  the  origin  and  radius  equal  to  ct/2.  For  a  multi-static 
radar,  the  bistatic  range  equation  has  to  be  used. 

tp  =  \I(X~XJ2  +(y-  y,xf  +  (Z  -  z,J2 1  c  +  yl(x-x„)2  +(y-  yj2  +  (z  -  zj2  (43) 

where  [xtx,  ytXt  ztx]  and  [xrx,  yrXi  zrx]  are  the  position  coordinates  of  the  transmitter  and 
receiver  respectively.  The  set  of  points  corresponding  to  the  same  reception  times  are 
ellipses  with  the  transmitter  and  receiver  locations  at  the  focal  points,  corresponding  to 
the  well-known  “ovals  of  Cassini,”  and  is  a  generalisation  of  the  case  of  monostatic  radar. 
For  a  moving  sensor,  the  exact  process  is  more  complicated  due  to  the  dynamic  change  of 
positions  from  range  sample  to  range  sample.  Ideally,  the  time  map  should  be 
recomputed  for  each  range  sample  before  being  used  to  update  the  image  map.  In  the 
implemented  model,  we  assume  that  the  time  map  does  not  change  significantly  during 
the  pulse  in  order  to  reduce  the  computational  load  and  simply  use  the  time  map  based  on 
the  sensor  positions  at  the  start  of  transmission.  This  assumption  is  likely  to  be  valid  as 
long  as  the  sensor  moves  much  less  than  the  required  imaging  resolution  cell  during  the 
pulse  duration. 
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The  implemented  model  makes  use  of  the  mapping  of  reception  time  to  each 
sample  point  in  the  desired  scene,  as  defined  by  a  set  of  contiguous  known  coordinates. 
This  map  of  reception  times  is  used  to  compare  against  reception  times  associated  with 
the  range  profiles  in  order  to  assign  these  range  profiles  values  to  the  right  spatial 
location.  Figure  16  illustrates  the  concept  of  this  algorithm,  which  is  described  as  follows: 

Step  1 :  Create  a  matrix  of  spatial  coordinate  points  representing  the  desired  scene. 
The  separation  between  each  of  the  points  shall  be  smaller  than  the  sampling  cell 
of  the  range  profile  data  to  be  used. 

Step  2:  Create  a  matrix  of  the  same  size  to  store  imaging  data  corresponding  to 
the  scene.  This  is  matrix  hereinafter  referred  to  as  the  “image  map.” 

Step  3:  Using  the  known  transmitter  and  receiver  location,  the  signal  reception 
time  after  travelling  to  each  spatial  coordinate  point  is  computed.  This  set  of 
reception  times,  each  corresponding  to  a  spatial  coordinate  in  the  desired  scene  is 
hereinafter  called  “time  map.” 

Step  4:  For  a  single  range  profile,  process  each  range  sample  (with  amplitude,  An 
and  reception  time,  trx,n)  as  follows: 

Step  4a:  Search  for  all  position  indices  in  the  time  map  with  receptions 
times  between  trxn  -7/2 and  trxn  +7J2. 

Step  4b:  Add  the  amplitude  value,  An,  to  the  corresponding  points  on  the 
image  map. 

When  completed,  the  set  of  { An }  in  the  image  map  forms  the  “image”  provided 
by  a  single  pulse. 

Step  5:  Repeat  Step  4  for  all  available  range  profiles  and  cumulatively  sum  the 
image  maps  from  all  the  range  profiles.  When  completed,  the  image  can  be 
visualized  by  plotting  the  magnitude  of  the  image  map. 

One  major  advantage  of  this  algorithm  is  that  it  automatically  performs  the 
polarimetric  (or  ellipsoidal)  interpolation  and  range  alignment  through  the  reception  time 
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mapping,  whereas  a  separate  step  is  typically  required  for  the  FFT  methods  to  work. 
Furthermore,  since  the  time  map  is  re-computed  based  on  sensor  configuration  for  each 
range  profile,  this  method  overcomes  the  limitation  present  in  both  the  FFT  and  Inverse 
Radon  Transform  methods  of  requiring  data  to  be  collected  at  uniformly  spaced  angles. 
This  method  is  also  very  flexible  in  the  sense  that  the  image  map  can  be  easily 
recomputed  for  a  different  scene  of  interest  using  the  same  range  profile  data  by  re¬ 
computing  the  time  map  corresponding  to  the  new  scene.  Furthermore,  different 
waveforms  may  be  used  in  the  generation  of  each  range  profile,  since  the  correlation 
receiver  takes  care  of  the  characteristics  of  each  waveform  when  generating  each  range 
profile. 


Image  Map 


Figure  16.  Illustration  of  the  backprojection  method 
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IV.  RESULTS 


In  this  chapter,  the  results  produced  by  the  implemented  model  are  presented. 
Source  codes  are  appended  in  Appendix  B  for  reference. 

A.  SIMULATION  PARAMETER  SELECTION 

The  key  limitation  of  the  simulation  model  is  its  heavy  computational  resource 
requirements,  in  terms  of  both  time  and  memory.  MATLAB  also  places  restrictions  on 
the  maximum  size  of  vectors  and  matrices,  which  limits  the  scenarios  that  can  be 
simulated.  In  this  Section,  we  describe  the  effect  of  simulation  parameter  choices  on  the 
simulations  so  that  appropriate  choices  may  be  made  regarding  scenarios. 

1.  Waveform  Definition 

The  waveform  parameters  to  be  defined  include  the  carrier  frequency,  pulse  width, 
pulse  bandwidth,  and  pulse  repetition  interval.  The  pulse  width  defines  the  time 
interval — according  to  Equation  (41) — for  which  samples  need  to  be  generated  for  each 
range  profile  of  a  point  target,  since  samples  outside  of  this  interval  will  be  zero.  The 
pulse  bandwidth  defines  the  necessary  sample  time  interval  to  recover  the  baseband 
signal  since  it  is  necessary  to  sample  the  pulse  sequence  minimally  at  the  Nyquist  rate  in 
order  to  avoid  aliasing.  In  order  to  ensure  that  the  sampling  interval  is  sufficient  taking 
into  account  effects  such  as  target  motion,  we  oversample  the  time  sequence  beyond 
Nyquist  sampling. 

The  number  of  samples  per  pulse  width  is  given  by  pulse  width  divided  by  the 
sampling  interval.  For  a  pulse  width  of  1  ps  and  a  pulse  bandwidth  of  150  MHz  (i.e.,  for 
1  meter  down-range  resolution),  a  sampling  rate  of  1  GHz  is  selected.  This  translates  to 
about  1000  samples  per  range  profile.  This  number  of  samples  is  multiplied  by  the 
number  of  pulses  to  be  simulated  per  run,  which  can  be  large  in  the  case  of  a  moving 
SAR.  Furthermore,  to  resample  the  received  pulse  to  obtain  evenly  spaced  samples,  the 
search  algorithm  is  processed  sample-by-sample  within  the  range  profile.  Thus,  the 
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simulation  of  a  long  pulse  and/or  a  pulse  with  high  bandwidth  will  lead  to  range  profiles 
with  large  numbers  of  time  samples,  and  hence,  require  long  simulation  times  and  also 
heavy  use  of  memory. 

2.  Imaged  Scene  Definition 

The  scene  to  be  imaged  in  the  simulation  is  defined  by  the  scene  center  location, 
down-range  width,  and  cross-range  width,  as  well  as  the  resolution  steps  in  down  range 
and  cross  range.  The  down  range  and  cross  range  widths,  together  with  the  respective 
resolution  steps,  defines  the  number  of  cells  in  the  image  map  used  in  the  backprojection 
algorithm.  The  resolution  steps  in  the  image  map  need  to  be  of  the  same  order  as  the 
range  profile  sampling  interval,  in  order  to  maintain  the  resolution  present  in  the  range 
profile  samples.  If  the  image  map  resolution  is  too  coarse,  the  range  profile  samples  may 
not  map  correctly  to  the  spatial  locations  in  the  image  map,  and  thus  affect  the  imaging 
outcome.  It  is  not  desirable  to  have  an  overly  fine  resolution,  as  this  greatly  increases  the 
demand  on  computational  resources.  For  example,  consider  the  100  x  100  meter  scene. 
With  a  sampling  rate  of  1  GHz,  the  range  sampling  interval  works  out  to  be  0.3  meters. 
This  means  that  the  image  map  is  a  334  x  334  matrix  (111556  samples).  In  the 
backprojection  algorithm,  the  updating  of  the  image  map  is  processed  on  a  sample -by¬ 
sample  basis  for  each  range  profile  to  obtain  the  image  map  for  a  single  pulse.  This 
process  is  repeated  for  all  range  profiles. 

We  conclude  this  section  by  noting  that,  while  the  implemented  algorithm  is 
sufficiently  flexible  to  allow  for  simulation  of  a  wide  range  of  scenarios,  judicious  choice 
of  simulation  parameters  is  necessary  to  manage  the  computational  load,  and  thus,  this 
consideration  limits  the  types  of  scenarios  that  can  be  simulated  using  this  model. 

B.  ONE-DIMENSION  MODEL  VERIFICATION 

This  section  focuses  on  the  verification  of  processes  leading  up  to  the  range 
profile  generation.  In  order  to  check  the  correctness  of  the  model,  the  outputs  from  the 
model  are  checked  against  known  theoretical  predictions.  The  results  from  these  checks 
are  consistent  with  theoretical  predictions.  This  lends  confidence  that  the  model  is 
accurate.  The  results  are  presented  as  follows: 
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1.  Doppler  Shift  for  Approaching  and  Receding  Target 

The  geometry  of  scenarios  involving  approaching  and  receding  targets  is 
illustrated  in  Figure  17.  The  theoretical  outcome  for  such  scenarios  is  well  known  from 
literature  and  was  also  presented  in  Chapter  EL  The  change  in  propagation  time  from 
sample  to  sample  should  correspond  to  the  predictions  from  Equation  (24),  when 
uniformly  spaced  transmission  samples  are  used  to  generate  the  reception  times.  Equation 
(24)  should  continue  to  be  valid  even  when  non-uniform  transmission  times  are  used  to 
obtain  the  required  uniformly  spaced  reception  times.  Figure  18  shows  the  results 
obtained  for  approaching  and  receding  targets  at  various  speeds.  As  expected,  the 
propagation  time  changes  as  a  linear  function  of  the  transmission  time,  with  the  slope 
determined  by  the  radial  velocity.  The  differences  obtained  between  the  obtained  Doppler 
(by  computing  the  gradient  of  the  propagation  time  change)  and  predicted  Doppler  shifts 
are  less  than  10-6  Hz. 
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Figure  17.  Geometry  of  the  approaching  target  (Top)  and  receding  target 

(Bottom)  scenarios 
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Change  in  propagation  time  versus  transmission  time 


Figure  18.  Change  in  propagation  time  as  a  function  of  transmission  time  for 
various  approaching  and  receding  radial  velocities. 


To  check  for  minute  non-linearities  that  cannot  be  seen  from  observing  the 
“straight  line”  graphs  observed  in  Figure  18,  the  difference  between  the  change  in 
propagation  time  and  a  linear  line  based  on  theoretical  predictions  is  computed  and 
plotted.  A  result  of  zero  across  all  transmission  times  should  be  obtained  if  the  results 
conform  to  theory.  If  the  relationship  were  non-linear,  then  a  non-linear  relationship 
should  be  apparent  from  the  curve.  The  results  from  this  “deviation  from  linear”  curve, 
shown  in  Figure  19,  look  fairly  noisy,  although  the  magnitude  of  the  deviation,  at 
1(T13  ps,  is  very  small. 
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Figure  19.  Difference  between  the  simulated  change  in  propagation  time 
versus  theoretically  predicted  linear  line  as  a  function  of  transmission  time 

Following  some  investigation,  it  was  confirmed  that  this  is  due  to  quantizing 
noise  from  double  precision  floating  point  arithmetic  used  in  MATLAB.  We  provide  a 
short  explanation  here.  MATLAB — and  in  general,  double  precision  floating  point 
calculations — have  a  computational  dynamic  range  of  around  15-17  digits.  This 
essentially  means  that  when  we  have  two  numbers  whose  difference  in  magnitude  is 
greater  than  1016,  the  smaller  number  cannot  be  represented  correctly.  This  is  a 
consequence  of  the  digital  representation  of  such  floating  point  numbers  in  64  bits. 
MATLAB’s  Symbolic  Computation  Toolbox  provides  a  function  that  implements 
Variable  Precision  Arithmetic  (VPA),  which  allows  for  the  computation  of  results  to 
beyond  double  precision  accuracy.  The  main  tradeoff  in  the  use  of  VPA  is  a  great 
increase  in  the  amount  of  computer  memory  and  simulation  times  required.  The  result  of 
a  short  study  on  the  use  of  VPA  versus  double  precision  arithmetic  is  provided  in 
Appendix  C  for  information.  These  accuracies  are  mainly  required  only  when  a 
combination  of  low  Doppler  speeds  and  very  short  sampling  intervals  are  involved.  As 
such,  we  have  scoped  our  simulation  scenario  to  avoid  long  simulation  times  due  to  the 
use  of  VPA.  The  use  of  VPA  is  not  further  considered  in  this  thesis. 
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2.  Non-linear  Relationship  Between  Transmission  and  Reception  Times 

The  main  cases  of  interest  for  study  using  the  model  are  moving  targets  which  are 
not  directly  ingressing  or  egressing.  Hence,  we  attempt  to  verify  the  model  with  a 
representative  scenario.  We  pick  the  case  of  a  single  crossing  target,  starting  at  [0;-20;0] 
km  and  ending  at  [0;+20;0]  km.  The  target  is  simulated  at  various  constant  speeds  along 
the  cross-range  axis.  The  radar  is  100km  from  the  origin.  Figure  20  illustrates  this 
scenario. 
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Figure  20.  Geometry  for  scenario  with  crossing  target.  This  is  for  the  case  of  a 

target  with  crossing  speed  of  7.8  km/s 

From  the  derivation  in  Chapter  II,  we  can  expect  the  initial  variation  in 
propagation  time  to  be  non-linear  a  function  of  transmission  time.  Figure  21  shows  the 
results  obtained,  which  are  as  expected.  We  can  observe  that  the  magnitude  of  the  change 
in  propagation  time  does  not  depend  on  the  speed,  but  rather  on  the  location  of  the  target 
at  the  point  of  interaction  as  expected.  As  a  check,  we  can  compute  the  difference  in 
slant  range  from  the  radar  between  the  start  position  and  closest  point  of  approach  and 
compute  the  difference  in  propagation  time.  This  difference  works  out  to  be  13.2ps, 
which  is  consistent  was  the  simulation  results. 
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Figure  21 .  Change  in  propagation  time  (referenced  to  first  transmitted  sample) 
for  crossing  targets  at  100  m/s  (Left)  and  1000  m/s  (Right). 

Figure  22  shows  the  change  in  propagation  time  over  the  period  of  a  500  ps  pulse 
transmission  period  for  a  100  m/s  target  crossing  at  closest  point  of  approach  100km 
away.  Based  on  these  results,  we  can  see  that  the  change  in  propagation  time  is  less  than 
a  femtosecond.  Compared  to  the  period  of  waves  at  microwave  frequencies,  this  variation 
could  be  insignificant,  and  thus  justifies  the  adoption  of  the  start-stop  approximation  for 
slow  moving  targets.  However,  when  a  longer  pulse  is  adopted,  the  non-linearity  will 
become  more  significant.  Figure  23  shows  the  change  in  propagation  time  over  the  period 
of  a  50  millisecond  pulse  for  a  250  m/s  crossing  at  closest  point  of  approach  50km  away. 
In  this  case,  the  deviation  is  a  hundred-th  of  a  nanosecond,  which  is  significant  compared 
to  the  period  of  a  10  GHz  wave.  Hence,  it  can  be  seen  that  for  longer  waveforms,  there 
may  be  a  further  need  to  examine  whether  the  start-stop  approximation  is  appropriate. 


Transmission  time  (millisec) 

Figure  22.  Change  in  propagation  time  (referenced  to  start  of  transmission 
sample)  for  a  500  psec  pulse  on  a  lOOm/s  crossing  target  at  100km  away. 
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Figure  23.  Change  in  propagation  time  (referenced  to  start  of  transmission 
sample)  for  a  50  millisec  pulse  on  a  250m/s  crossing  target  at  50km  away. 


3.  Ambiguity  Functions  for  Known  Waveforms 

The  ambiguity  functions  of  waveforms  are  useful  tools  for  analyzing  the 
characteristics  of  a  radar  waveform  in  the  spatial  and  Doppler  domains.  To  check  the 
reception  process,  we  generated  the  ambiguity  functions  for  a  monostatic  setup  by 
allowing  the  received  signal  to  be  the  same  as  the  transmitted  signal  and  passing  this 
thought  the  matched  filter.  This  is  done  for  a  rectangular  pulse,  a  Barker- 13  binary  phase 
coded  pulse,  and  a  LFM  pulse.  Figure  24,  Figure  25,  and  Figure  26  show  the  ambiguity 
functions  obtained.  These  compare  well  with  known  ambiguity  function  plots  found  in 
literature  ([27],  pp.55,  pp.  1 14  and  pp.58). 


Figure  24.  Ambiguity  function  and  associated  zero  Doppler  cut  for  a 
rectangular  pulse  from  model.  The  delay  and  frequency  axes  have  been 
normalized  using  the  simulated  pulse  width  and  pulse  bandwidth  (=  reciprocal  of 
pulse  width  in  this  case),  respectively. 
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Figure  25.  Ambiguity  function  and  associated  zero  Doppler  cut  for  a  Barker 
13  binary  phased  coded  pulse  from  model.  The  delay  and  frequency  axes  have 
been  normalized  using  the  subpulse  width  and  pulse  bandwidth,  respectively. 


Figure  26.  Ambiguity  function  and  associated  zero  Doppler  cut  generated  for 
a  LFM  pulse  (Bandwidth=15MHz)  from  model.  The  delay  and  frequency  axes 
have  been  normalized  using  the  pulse  width  and  pulse  bandwidth,  respectively. 


C.  TWO-DIMENSION  IMAGING  MODEL  VERIFICATION 

This  section  focuses  on  the  verification  of  processes  related  to  image 
reconstruction. 

1.  Time  Map  Implementation  for  2-D  Imaging 

The  correct  generation  of  the  time  map  is  critical  for  image  reconstruction  using 
the  backprojection  method.  Figure  27  and  Figure  28  show  the  time  map  generated  for 
both  monostatic  and  bistatic  radar  configurations.  These  maps  were  generated  with  coarse 
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spatial  sampling  at  1km  intervals  to  provide  a  “big  picture”  view  of  the  time  map.  As 
expected,  the  monostatic  radar  presented  a  circular  time  map  and  the  multistatic  radar 
provided  an  ellipsoidal  map. 


Figure  27.  Contour  plot  of  time  map  for  a  monostatic  radar  located  at 

[20km, 20km] . 


Figure  28.  Contour  plot  of  time  map  for  bistatic  radar  with  transmitter  at 

[50km, -50km]  and  receiver  at  [10km, -10km]. 
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Figure  29  shows  the  typical  time  map  that  will  be  used  in  the  imaging  process. 
This  is  a  detailed  time  map  generated  for  a  scene  of  400m  x  400m  in  size  at  10cm 
intervals.  One  observation  here  is  that,  although  the  time  map  looks  fairly  elliptical  in 
Figure  28,  the  contour  plots  here  look  fairly  linear  even  though  the  scene  center  is  only 
about  14km  away  from  the  receiver.  This  means  that  there  could  be  potential  for  linear 
mapping  between  the  spatial  and  time  domains  over  such  small  scenes  that  could  simplify 
the  implementation  of  backprojection  algorithms. 


Figure  29.  Detailed  time  map  for  400m  x  400m  area  around  scene  center 
based  on  the  geometry  shown  in  Figure  28.  Contours  are  set  at  5ps 
intervals  for  display  only.  Samples  are  generated  at  10cm  intervals. 

2.  Single  Pulse  Image  Map 

The  image  map  due  to  a  single  pulse  is  checked  to  ensure  proper  mapping 
between  the  time  map  results  and  the  image  map.  A  target  is  placed  at  [-50m, 50m]  for 
this  check.  Figure  30  and  Figure  31  show  the  image  map  generated  using  the  monostatic 
and  bistatic  sensor,  respectively. 
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Figure  30.  Single  LFM  pulse  (15  MHz  bandwidth)  image  map  using 
monostatic  configuration,  with  target  at  [-50m, 50m]  using  the  backprojection 

algorithm. 
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Figure  31.  Single  LFM  pulse  (15  MHz  bandwidth)  image  map  using  bistatic 
configuration,  with  target  at  [-50m, 50m]  using  the  backprojection  algorithm. 
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The  combined  image  map  is  shown  in  Figure  32.  It  can  be  seen  that  the  target  is 
clearly  localized  at  its  true  location  with  just  two  pulses  in  this  case. 


4\. 


Figure  32.  Combined  image  map  formed  from  the  outputs  of  the  monostatic 
and  bistatic  radars.  The  target  is  clearly  localized  at  its  true  location,  [-50m,  50m]. 

D.  EFFECT  OF  TARGET  SPEED  ON  RANGE  PROFILE 

Cheney  and  Borden  [15]  listed,  as  a  guideline,  a  validity  domain  for  start  stop 
approximation  as  v  <  X  /  Tw  where  v  is  the  relative  speed  between  sensor  and  target,  X  is 
the  wavelength  at  the  center  frequency  and  Tw  is  the  pulse  width.  This  essentially  means 
that  the  motion  between  sensor  and  target  over  a  single  pulse  width  is  less  than  a 
wavelength  and  sets  limits  on  the  desired  target  speed  or  pulse-width,  given  that  the 
transmission  frequency  is  typically  chosen  based  on  other  considerations  previously 
discussed.  The  effects  within  the  validity  domain  are  examined  here. 

As  a  reference,  a  single  LFM  pulse  with  pulse  width  of  lps  and  bandwidth  of  15 
MHz  is  used.  A  carrier  frequency  of  1  GHz  is  chosen,  thus  the  wavelength  is  30  cm.  The 
radar  is  100km  from  the  target.  The  effects  on  the  range  profile  are  observed  when  the 
speed  for  an  inbound  target  is  varied  from  lkm/s  to  lOOkm/s,  or  equivalently  v  varies 
from  0.003  X  /  Tw  to  0.3  X  /  Tw .  The  range  profile  is  selected  at  the  zero  Doppler  cut  of  the 
correlation  filter  bank.  Figure  33  shows  the  results. 
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Figure  33.  Effect  of  target  speed  (inbound  target)  from  1  km/s  to  100  km/s  on 
the  range  profile  output  for  LFM  pulse  with  Tw  =  lps,  bandwidth  =15  MHz.  The 
delay  axis  is  normalized  by  the  pulse  width,  the  amplitude  is  normalized  using  the 
peak  of  the  auto-correlation  function. 

It  can  be  clearly  seen  that  the  range  profile  peak  is  shifted  from  the  “correct” 
range  in  the  profile  and  the  peak  amplitude  is  reduced  by  more  than  3dB  even  when  the 
target  speed  is  still  much  smaller  than  A  /  Tw .  These  effects  will  contribute  to  artifacts  in 
the  radar  imaging  processes. 

Figure  34  shows  the  output  for  the  same  waveform  with  a  lOOkm/s  target  when 
the  radar  is  at  45  degree  aspect  angle  to  the  target  and  a  target  with  velocity  at  90  deg  to 
the  down-range  direction.  It  can  be  seen  that  the  range  profile  output  is  less  impacted 
even  though  the  speed  is  the  same.  This  is  because  the  target  is  not  approaching  head-on, 
and  so  the  radar  waveform  is  experiencing  correspondingly  less  Doppler  modulation. 
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Figure  34.  Effect  on  the  range  profile  output  for  LFM  pulse  with  Tw  =  1  ps, 
bandwidth  =15  MHz  when  a  lOOkm/s  target  is  at  45  degree  aspect  angle  (Left) 
and  when  velocity  is  at  90  deg  to  the  down-range  direction  (Right).  The  delay  axis 
is  normalized  by  the  pulse  width,  the  amplitude  is  normalized  using  the  peak  of 

the  auto-correlation  function. 


The  effects  of  changing  the  LFM  pulse  are  examined.  Figure  35  shows  the  output 
with  an  inbound  target  at  lOOkm/s,  with  a  longer  pulse  width  of  Tw  =  lOps,  bandwidth  = 
15  MHz.  It  can  be  seen  that,  in  both  cases,  the  distortions  introduced  by  the  target’s 
motion  are  reduced.  In  the  case  of  a  longer  pulse,  it  can  be  seen  that  the  actual  range  shift 
introduced  by  the  Doppler  shift  remains  about  the  same  as  the  reference,  since  the  pulse 
width  has  increased  ten-fold. 


Figure  35.  Effect  on  the  range  profile  output  for  LFM  pulse  with  longer 
pulse, Tw  =  lOps  (Left)  and  higher  bandwidth  =150  MHz  (Right)  with  an  lOOkm/s 
inbound  target.  The  delay  axis  is  normalized  by  the  pulse  width,  the  amplitude  is 
normalized  using  the  peak  of  the  auto-correlation  function. 
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To  explore  the  dependence  on  waveforms,  the  effects  of  the  same  scenario  on  a 
Barker- 13  binary  phase  coded  signal  are  examined.  The  Barker  code  is  often  stated  in 
literature  as  having  poor  “Doppler  tolerance.”  [27]  The  bandwidth  of  this  waveform  is 
kept  at  15  MHz  to  maintain  the  same  range  resolution;  this  gives  a  slightly  shorter  pulse 
width  of  0.867  ps.  Figure  36  shows  the  effect  of  target  motion  on  this  waveform  at  higher 
target  speeds  of  50km/s  and  lOOkm/s.  It  can  be  seen  that  the  effect  of  target  motion 
essentially  is  to  raise  the  range  sidelobes  and  lower  the  peak  response  of  the  return,  but 
does  not  introduce  a  range  error.  This  is  likely  to  cause  artifacts  that  are  different  from 
that  of  LFM  pulses  in  the  imaging  outcome — for  example,  the  raised  sidelobes  from  large 
targets  masking  out  other  nearby  features  in  the  image. 


Figure  36.  Effect  on  the  range  profile  output  for  Barker-13  binary  phase 
coded  pulse  with  subpulse  width,  Tw  =  0.667 ps  for  a  50  km/s  and  100  km/s 
inbound  target.  The  delay  axis  is  normalized  by  the  subpulse  width,  the  amplitude 
is  normalized  using  the  peak  of  the  auto-correlation  function. 

Hence,  it  can  be  seen  that  the  artifacts  due  to  imaging  will  be  dependent  on  the 
characteristics  of  the  waveform  chosen,  in  addition  to  the  target  characteristics. 

E.  MULTI-STATIC  RADAR  CASE  STUDY 

In  this  section,  we  present  the  imaging  outcomes  for  two  different  configurations 
of  multistatic  radar  systems.  The  first  is  a  collection  of  monostatic  radars,  assumed  to  be 
perfectly  synchronized;  and  the  second,  a  true  multi-static  system  where  the  transmitter 
and  receiver  are  not  collocated. 
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1.  Multiple  Monostatic  Radars 

This  configuration  consists  of  four  monostatic  radars  arranged  as  follows:  Radar  1 
is  at  [-100,100,0]  km,  Radar  2  is  at  [0,100,0]  km,  Radar  3  is  at  [50,0,0]  km,  and  Radar  4 
is  at  [50,86.6,0]  km.  This  corresponds  to  observation  angles  of  -30,  0,  45,  and  90  degrees, 
but  at  differing  ranges.  LFM  pulses  with  pulsewidths  of  lps  and  bandwidth  of  15  MHz 
are  used  in  all  the  radars. 

A  point  target  located  at  initial  position  [-50,50,0]  meters  is  observed.  The  sensor 
and  target  layout  is  shown  in  Figure  37.  The  imaging  outcome,  i.e.,  Point  Spread 
Function  (PSF),  for  a  stationary  target  is  shown  in  Figure  38.  The  target  is  well  localized 
at  its  true  position  with  the  expected  resolution  from  a  15  MHz  bandwidth  (i.e.,  10 
meters).  The  lines  of  backprojection  from  the  radar  positions  are  clearly  observable. 
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Figure  37.  Layout  of  radars  for  multiple  monostatic  imaging  case  study 
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Figure  38.  Point  Spread  Function  for  stationary  target  at  [-50,50]  meters. 

Radars  are  at  [-100, 100,0] km,  [0,100,0]km,  [50,0,0]km  and  [50,86.6,0]km. 

Amplitudes  are  normalized  by  the  number  of  pulses  processed. 

The  effect  of  target  motion  is  now  investigated.  Figure  39  shows  the  PSF  for  a 
target  moving  at  7.5  km/s  parallel  to  the  Y  and  X-axes.  A  close-up  view  of  the  imaging 
PSF  for  the  stationary  target  is  provided  for  reference.  We  note  that  during  the  interaction 
with  the  pulse,  the  target  moves  less  than  1cm.  It  can  be  seen  that  the  imaging  responses 
are  clearly  shifted  away  from  the  true  target  location  by  approximate  2.5  meters  in  each 
case.  Visually,  it  can  be  seen  that  the  shape  of  the  PSF  has  broadened  slightly  across  the 
axis  of  motion.  The  peak  responses  in  both  cases  have  also  decreased  compared  to  the 
stationary  target. 

A  hypothetical  high  speed  target  travelling  at  50km/s  in  the  direction  of  the  +Y 
axis  has  been  used  to  illustrate  its  impact  on  the  imaging  PSF.  The  target  moves  5cm 
during  the  lps  pulse  duration.  The  observed  PSF  is  shown  in  Figure  40.  This  image 
clearly  shows  the  same  degradations  shown  in  the  previous  figures,  including  a  wrong 
target  location  and  decrease  in  the  peak  responses.  In  addition,  it  can  be  clearly  observed 
that  the  backprojected  range  profiles  fail  to  converge  at  the  same  focal  point,  which 
ideally  should  be  at  the  true  location  of  the  target.  In  effect,  the  target  has  become 
defocused,  resulting  in  multiple  spurious  peak  responses.  This  is  because  different  radars 
“see”  different  radial  velocities  due  to  their  aspect  angle,  and  thus  the  range  error 
introduced  varies  between  different  range  profiles.  When  this  deviation  is  sufficiently 
large,  the  projection  lines  will  no  longer  intersect  at  the  same  location. 
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Figure  39.  Point  Spread  Function  for  target  with  initial  position  [-50,50,0] 
meters  moving  at  7.5  km/s,  (a)  parallel  to  +Y-axis  (Top  Left),  (b)  parallel  to  +X- 
axis  (Top  Right).  The  PSF  for  the  stationary  target  (Bottom)  is  shown  for 

reference. 
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Figure  40.  Point  Spread  Function  for  target  with  initial  position  [-50,50,0] 

meters  moving  at  50  km/s. 
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The  same  simulations  are  repeated  on  a  set  of  radars  using  Barker-13  binary  phase 
coded  waveforms.  As  expected,  it  can  be  seen  that  the  range  sidelobes  appear  more 
prominently  in  the  imaging  results.  The  results  are  shown  in  Figure  41.  The  observed 
effects  are  similar  to  the  case  of  the  LFM  waveform. 

2.  Multistatic  Radars 

In  this  case,  we  check  an  actual  multi-static  configuration  where  the  transmitter 
and  receiver  are  not  co-located.  The  locations  of  the  receivers  are  the  same  as  the  radar 
locations  in  the  previous  scenario:  [-100,100,0]km,  [0,100,0]km,  [50,0,0]km,  and 
[50,86.6,0]  km.  The  transmitter  is  located  at  [0,100,0]km  and  transmits  LFM  pulses  with 
pulsewidths  of  lps  and  bandwidth  of  15  MHz.  Again,  the  target  is  placed  at  [-50,50,0] 
meters.  Figure  42  shows  the  layout  of  the  scenario  and  the  PSF  for  a  stationary  target. 
Again,  the  target  is  well  localized  at  its  true  position  with  the  expected  resolution  from  a 
15  MHz  bandwidth  (i.e.,  10  meters).  The  lines  of  backprojection  from  the  radar  positions 
are  clearly  observable  and  clearly  different  in  orientation  from  that  of  the  multiple 
monostatic  radar  case.  This  is  because  the  backprojections  of  the  multistatic  radar 
wavefronts  are  aligned  to  ellipses  with  the  transmitter  and  receiver  as  the  foci,  unlike 
circular  arcs  in  the  monostatic  radar  case. 
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Figure  41.  Point  Spread  Function  for  configuration  using  Barker- 13  binary 
phase  coded  signal  for  target  towards  +Y  direction  at,  (a)  7.5  km/s  (Top  Left),  (b) 
50  km/s.  The  point  spread  function  for  a  stationary  target  (Bottom)  is  shown  for 

reference. 
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Figure  42.  Layout  of  target,  transmitter  and  receivers  (Left)  and  Point  Spread 
Function  (Right)  for  the  multistatic  radar  scenario. 
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It  can  be  seen  from  Figure  43  that  the  effects  of  a  moving  target  are  similar  to  the 
multiple  monostatic  radar  case. 
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Figure  43.  Multistatic  configuration  Point  Spread  Function  for  target 
with  initial  position  [-50,50,0]  meters  moving  at  7.5  km/s,  (a)  parallel  to 
+Y-axis  (Top  Left),  (b)  parallel  to  +X-axis  (Top  Right).  The  PSF  for  the 
stationary  target  (Bottom)  is  shown  for  reference. 


We  conclude  this  section  by  summarizing  the  observations  that  have  been  made. 
It  has  been  observed  that  moving  targets  may  introduce  a  shift  in  the  location  of  the  target 
image.  This  follows  from  the  distortions  in  the  range  profile  that  have  been  observed  in 
the  previous  section.  When  the  target  is  sufficiently  fast,  the  range  shift  introduced  varies 
between  different  range  profile.  This  deviation  may  cause  standard  backprojection  to  fail, 
since  the  backprojected  range  profiles  will  no  longer  intersect  at  the  same  point,  leading 
to  image  defocusing. 
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F.  MONO-STATIC  SINGLE  CHANNEL  SAR  CASE  STUDY 

In  this  section,  we  examine  the  effects  of  a  moving  target  on  single  channel  SAR. 

1.  SAR  Cross  Range  Resolution 

Prior  to  observing  the  effects  of  a  moving  target,  the  SAR  model  is  checked  to 
verify  that  the  model  conforms  to  theoretical  predictions  for  known  cases.  A  LFM 
waveform  with  pulse-width  of  lps  and  modulation  bandwidth  of  15  MHz  (corresponding 
to  down-range  resolution  of  10  meters)  is  used  in  this  series  of  verifications.  A  carrier 
frequency  of  100  MHz,  with  sampling  frequency  of  1  GHz,  is  chosen  for  a  manageable 
computation  load.  The  expected  cross-range  resolution  of  this  model  can  be  computed 
using  Equation  (4).  When  processing  the  SAR  imagery,  a  salient  point  to  note  is  that  the 
SAR  data  needs  to  be  processed  at  the  carrier  frequency  in  order  to  achieve  the  full 
resolution  as  predicted  [13]. 

The  effect  of  observation  width,  L,  (i.e.,  SAR  speed  multiplied  by  observation 
time)  is  first  examined.  As  a  reference  case,  L  is  chosen  to  be  3km.  The  SAR  is 
positioned  with  a  down-range  of  10km,  with  a  speed  of  50  m/s  (i.e.,  100  kts).  The  SAR 
path  is  chosen  such  that  the  target,  located  at  the  origin,  appears  broadside  in  the  middle 
of  the  observation  width  and  samples  are  taken  at  half-wavelength  intervals  in  the 
synthetic  aperture  domain.  This  geometry  is  illustrated  in  Figure  44.  With  these 
parameters,  the  PSF  obtained  is  shown  in  Figure  45.  It  can  be  seen  that  the  cross-range 
resolution  achieved  is  around  10  meters  (null-to-null),  which  is  as  expected. 

The  PSF  is  generated  for  observation  widths  of  1.5km  and  6km.  The  results 
shown  correspond  well  to  theory,  with  cross  range  resolutions  (null-to-null)  of  20m  and 
5m,  respectively.  The  PSF  is  also  observed  for  carrier  frequencies  of  75MHz  and 
150MHz,  using  the  reference  observation  width  of  3km.  The  results  shown  in  Figure  47 
again  correspond  well  to  theory  with  cross  range  resolutions  (null-to-null)  of  13.3  meters 
and  6.7  meters,  respectively.  These  verifications  show  that  the  implemented  model  is 
consistent  with  theory,  in  terms  of  both  the  expected  down-range  and  cross-range 
resolutions. 
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Figure  44.  Geometry  of  reference  SAR  scenario.  The  SAR’s  speed  is  chosen 
to  be  50m/s.  The  X  and  Y  axes  are  the  down-range  and  cross-range  directions, 

respectively. 
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Figure  45.  Point  Spread  Function  for  reference  SAR  scenario. 


Down  Range  (meters)  Down  Range  (meters) 


Figure  46.  Point  Spread  Function  for  SAR  scenario  with  observation  widths 

of  1.5km  (Left)  and  6km  (Right) 
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Figure  47.  Point  Spread  Function  for  SAR  scenario  with  carrier  frequencies  of 

50  MHz  (Left)  and  150  MHz  (Right) 


2.  Point  Spread  Function  for  a  Moving  Target 

Here,  the  PSF  is  observed  for  a  moving  target.  Figure  48  and  Figure  49  show  the 
PSFs  for  a  SAR  scenario  with  target,  initially  located  at  [0,0,0]  meters,  moving  at  0.1  m/s 
in  the  downrange  (positive  X)  and  cross-range  (positive  Y)  directions,  respectively.  The 
targets  move  0.1pm  during  a  single  pulse  duration  and  6  meters  during  the  entire 
observation  period.  This  slow  speed  is  chosen  to  illustrate  cases  where  there  is  very  little 
Doppler  modulation,  which  results  in  negligible  range  shift  of  the  range  profile’s  peak 
output  compared  to  a  stationary  target.  It  can  be  clearly  seen  in  both  cases  of  down-range 
and  cross-range  motions,  that  there  are  errors  introduced  in  the  location  of  the  target  and 
the  peak  output  of  the  responses  have  slightly  decreased.  For  the  target  moving  in  the 
downrange  direction,  the  error  in  the  cross-range  location  is  20  meters  as  predicted  with 
the  formula  in  [16]. 

In  the  case  of  the  cross-range  motion,  the  resolution  of  the  PSF  has  been  degraded, 
as  the  original  main  lobe  of  the  PSF  has  merged  with  its  adjacent  sidelobes  to  create  a 
broadened  response.  This  is  the  defocusing  effect  that  is  discussed  in  [16]  and  [17]. 
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Figure  48.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  0.1  m/s  in  the  down-range  direction.  The  PSF 
for  a  stationary  target  is  shown  (Bottom)  for  reference. 


Figure  49.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  0.1  m/s  in  the  cross-range  direction. 


We  now  examine  the  PSF  for  a  faster  target  which  is  moving  at  lOm/s.  Figure  50 
shows  the  PSF  for  a  SAR  scenario  with  a  target  initially  located  at  [0,0,0]  meters,  moving 
at  10  m/s  in  the  downrange  (positive  X)  and  cross-range  (positive  Y)  directions.  The 
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target  moves  1 00 jam  during  a  single  pulse  duration  and  600  meters  during  the  entire 
observation  period.  This  speed  is  chosen  to  illustrate  cases  where  the  target  moves  a 
distance  which  is  much  greater  than  the  imaging  resolution  during  the  observation  period. 
The  peak  response  is  again  shifted  in  cross-range  by  the  amount  predicted  using  [16].  We 
see  that  the  PSF  is  smeared  out,  the  peak  response  is  greatly  reduced,  and  there  are 
multiple  peaks  in  the  response.  This  is  due  to  the  fact  that  the  target  continuously  varies 
its  physical  location  throughout  the  observation  period,  resulting  in  an  effect  similar  to 
“motion  blur”  in  optical  cameras. 
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Figure  50.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  10  m/s  in  the  down-range  direction. 


Figure  51.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  10  m/s  in  the  cross-range  direction. 
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The  PSF  for  targets  moving  at  lOm/s  is  now  further  examined  with  a  faster  SAR 
speed  of  250m/s  to  reduce  the  observation  time.  Target  motion  during  a  single  pulse 
duration  remains  the  same,  but  the  target  now  moves  only  120  meters  during  the  entire 
observation  period.  Figure  52  and  Figure  53  show  the  PSF  for  a  SAR  scenario  with  target, 
initially  located  at  [0,0,0]  meters,  moving  at  10  m/s  in  the  downrange  (positive  X)  and 
cross-range  (positive  Y)  directions,  respectively.  It  can  be  observed  that  the  smearing  and 
cross-range  error  is  greatly  reduced  in  both  cases  despite  the  greater  relative  motion 
between  the  SAR  and  the  target.  This  is  consistent  with  Raney’s  [16]  conclusion  that  a 
large  SAR  speed  reduces  the  artifacts  induced  by  target  motion. 


Cross  Range  (meters)  Do*"  Ran9e  (meters) 


Figure  52.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  10  m/s  in  the  down-range  direction  and  with 

SAR  moving  at  250  m/s. 
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Figure  53.  Point  Spread  Function  for  SAR  scenario  with  target,  initially 
located  at  [0,0,0]  meters,  moving  at  10  m/s  in  the  cross-range  direction  and  with 

SAR  moving  at  250  m/s. 
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In  summary,  the  effects  of  target  motion  on  SAR  imaging  have  been  examined. 
The  phase  distortions  introduced  by  a  motion  target  will  lead  to  an  error  in  the  cross¬ 
range  location  of  the  target  and  defocusing.  In  addition,  when  the  target  moves  over  a 
long  distance  during  the  observation  period,  the  resultant  imaging  response  is  smeared.  In 
such  cases,  a  higher  SAR  speed  will  help  to  reduce  the  artifacts  caused  by  target  motion. 
However,  in  practical  systems,  higher  SAR  speeds  will  require  the  use  of  higher  Pulse 
Repetition  Frequencies  (PRF)  to  maintain  the  sampling  distance  in  the  synthetic  aperture 
domain.  There  are  other  practical  constraints  on  PRF  selection  considerations,  such  as  the 
required  radar  unambiguous  range,  as  well  as  radar  hardware  limitations.  These  other 
considerations  have  not  been  factored  into  our  simulations,  but  will  nevertheless  need  to 
be  addressed  in  a  real  system. 
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V.  CONCLUSIONS 


In  this  thesis,  we  study  the  imaging  degradations  or  artifacts  that  occur  when  a 
scene  to  be  imaged  by  radar  contains  moving  targets. 

The  physics  of  interaction  between  radar  waves  and  moving  targets  were  studied, 
and  a  simulation  was  developed,  using  MATLAB  to  model  the  resultant  signal  received. 
A  backprojection  algorithm  was  also  implemented  for  radar  image  reconstruction  from 
multiple  one-dimensional  range  profiles,  each  corresponding  to  the  output  from  a  single 
radar  pulse.  The  model  implemented  is  flexible  and  provides  the  means  to  study  the 
effects  of  target  motion  at  different  stages  in  the  image  formation  process,  and  also  for 
different  radar  configurations,  including  multistatic  radars  and  Synthetic  Aperture  Radar 
(SAR).  Most  importantly,  this  model  also  does  not  rely  on  the  commonly  used  start-stop 
approximation  which  assumes  that  the  target  and  sensor  are  both  stationary  throughout 
the  transmission  and  reception  process. 

Validations  were  performed,  with  results  checked  against  established  theories  or 
literature  to  provide  confidence  in  the  model.  The  key  limitation  of  this  model  is  its  high 
computational  resource  requirements  for  high  bandwidth  or  long  pulses,  due  to  the 
sample  sequence  lengths  that  need  to  be  simulated  for  proper  results. 

Through  these  studies,  it  was  noted  that  there  are  a  few  main  effects  that  lead  to 
the  degradation  of  a  radar  image.  The  first  effect  is  the  distortion  introduced  to  range 
profiles  due  to  the  mismatch  of  the  received  waveform  vis-a-vis  the  “matched  filter.” 
This  exact  effect  is  waveform-dependent,  but  in  the  commonly  used  Linear  Frequency 
Modulated  (LFM)  waveforms,  other  than  a  broadening  of  the  main  lobe  which  results  in 
smearing  in  the  final  image,  a  range  shift  is  introduced.  This  leads  to  the  object  being 
placed  in  the  wrong  location  on  the  final  image,  commonly  known  as  the  “train  off  the 
tracks”  or  “ship  off  its  wake”  phenomenon  by  SAR  image  analysts.  The  second  effect, 
evident  from  the  studies  on  multi-static  imaging,  stems  from  the  range  error  from  range 
profiles  affecting  the  intersection  points  of  the  back-projection  algorithm,  resulting  in 
defocusing.  Lastly,  for  a  moving  SAR  platform,  as  the  samples  are  taken  from  different 
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instants  in  time,  the  actual  physical  position  of  the  target  varies  throughout  the  imaging 
process.  In  the  case  of  slow  targets,  this  leads  to  defocusing  and  positional  errors, 
especially  in  the  cross-range  domain.  For  faster  targets,  the  intersections  of  the 
backprojections  do  not  intersect  correctly,  even  if  the  location  error  due  to  the  Doppler 
shift  has  been  corrected,  resulting  in  a  smeared  imaging  response.  This  is  similar  to 
“motion  blur”  experienced  in  optical  cameras  with  a  fast  object.  A  faster  SAR  speed  may 
help  to  reduce  these  artifacts.  These  observed  effects  are  consistent  with  descriptions 
from  literature,  but  have  been  observed  in  detail  due  to  the  flexibility  of  the  implemented 
model. 


The  implemented  model  is  sufficiently  flexible  to  allow  for  the  study  of  different 
waveforms  distorted  by  propagation,  without  first  having  to  derive  closed-form 
expressions,  which  are  typically  difficult  to  derive  in  the  moving  target  and  multi-static 
radar  cases.  Hence,  the  model  is  useful  for  further  studies  into  artifacts  introduced  by 
motion  in  radar  imaging,  as  well  as  methods  to  address  these  artifacts.  Possible 
opportunities  for  future  work  to  extend  or  apply  this  model  include  the  following: 

•  Using  the  model  to  further  study  aspects  of  radar  imaging  under  motion 
using  different  waveforms  and/or  matched  filter  generation. 

•  Using  the  model  to  generate  matched  filters  across  both  spatial  and 
velocity  domains  corresponding  to  the  Borden  Cheney  procedure  so  that 
its  performance  on  SAR  data  can  be  analyzed. 

•  Using  the  model  to  investigate  the  effects  of  time  synchronization  errors 
on  the  imaging  outcomes  of  multistatic  radars.  This  could  possibly  be 
done  by  adding  such  errors  to  the  time  map  for  image  formation. 

•  Developing  a  model  to  exactly  represent  the  “time  map”  of  moving  targets 
and  study  the  “time-map”  distortions  introduced  due  to  motion  vis-a-vis 
static  sensors  and  targets. 

•  Attempt  the  derivation  of  similar  expressions  as  Equations  (13)  and  (14) 
for  circular  flight  paths  and  geo-centric  circular  orbits  (e.g.,  useful  for 
satellites) 
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APPENDIX  A.  DERIVATION  OF  EQUATIONS 


This  appendix  provides  the  detailed  derivations  of  key  equations  that  were 
presented  in  Chapter  II. 

A.  PROPAGATION  TIME  FROM  TRANSMITTER  TO  TARGET 

Definition  of  terms: 

A tjX  =  one  way  propagation  time  from  transmitter  to  target 
t,x  =  transmission  time 

Target  position  at  t  =  0:  utgl  =  [uxtgt ;  uytgt ;  uz  tgt  ] 

Target  velocity  (constant):  vtgt  =b’XJgfWyJtgtWlJtgt\ 

Transmitter  position  at  t  =  0:  iitx  =  [nx  a ;  u  tx ;  u _  n  ] 

Target  velocity  (constant):  Vlx  =[vxa;vya;vztx] 


Since  the  wave  has  the  same  speed  in  all  directions,  the  wave  intercepts  the  target  when 
their  ranges  are  equal  from  the  transmitter.  From  the  frame  of  the  transmitter,  we  have 

REMwaVe  (A  +  )  =  cK  and  Rtgt  (ttt  +  A tix )  =  \ utgt  (ttx )  +  vtglA tix  -  utx  (tlx  )| 


Equating  the  two  ranges: 


Rtgt  i}tx  +  ^ix  )  R EM  wave  (C  +  ^ix  ) 

=>  cAtix  =  (A ) + ~  K  (A ) 
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=>  C"A4  =  \U,gt  (4  )  +  V'grAtu  ~  U,x  (4  )f 


Define  ut  (ta )  =  utgt  (0)  +  vtgtttx  -  utx  (0)  -  vtg,ttx 
=> c2Atfx  =  (u,  (ttx )  +  vtglA 4 ) •  (n,  (ttx)  +  vtgtA 4 ) 

=>  0  =  \«,  (ttx  )|2  +  2at  (ttx  )  •  v,gtAtix  +  (|v*  r-c2)  A  e 


Solving  the  quadratic  equation  for  tp  and  choosing  the  causal  solution  where  Atix  >  0 


K(ttx)  = 


-2Ut(ta)-vlg,± 


V 


4(4(4)-^/)  -4 [\y. 


2  ./ 2 


-c1  \u,(tj\2 


21 IV  -4 


4(4) ’V  +< 


(  »,  (4)  '  V'gt  )2  +  (c2  - \vv  f  )  |4  (4) | 


C2  -IV 


tgt 


This  yields  Equation  (13)  from  Chapter  II. 


B.  PROPAGATION  TIME  FROM  TARGET  TO  RECEIVER 

Definition  of  terms: 

Atrx  =  one  way  propagation  time  from  target  to  receiver 
Atix  =  one  way  propagation  time  from  transmitter  to  target 
ta  =  transmission  time 

Target  position  at  t  =  0:  utgt  =[uXJgt;uyJgt\utJgt] 

Target  velocity  (constant):  vtgt  =  [vxtgt ;  VyW ;  Vztgt  ] 

Transmitter  position  at  t  =  0:  Utx  =  [ux  tx ;  u  a ;  Uztx  ] 

Target  velocity  (constant):  vtx  =  [v xtx\V  ytx\V ztx\ 

Receiver  position  at  t  =  0:  urx  =  [uxrx\Uyrx\Uzrx] 

Receiver  velocity  (constant):  V rx  =  [vx  rx ;  v  rx ;  V,  rx  ] 
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The  key  difference  from  the  case  in  Appendix  A.I  is  that  the  target  is  now  the  source  of 
the  EM  wave,  which  it  emits  at  time  ttx+Atix,  and  the  receiver  is  the  object  to  be 
intercepted  by  the  wave,  which  is  received  at  trx  =  ttx+Atix+Atrx.  Since  the  wave  has  the 
same  speed  in  all  directions,  the  wave  intercepts  the  target  when  their  ranges  are  equal 
from  the  target. 


From  the  frame  of  the  target,  we  have 

Rm.  ( t  )  =  cAt  and  R  At  )  =  |m  (t.  +A  t.)  +  v  At  -ii.At.  +  At  )\ 

EMwave  \  rx )  rx  tgt  \  rx  J  \  rx  \  tx  ix  J  rx  rx  tgt  \  tx  ix )  | 


Equating  the  two  ranges: 


K  (t  r\  )  R EMwave  A rx) 

=>  cAt  —  u  (t,  +  At.  )+v  At  —u,  At,  +  At.  )| 

rx  rx  \  tx  ix  J  rx  rx  tgt  \  tx  ix  J  \ 

i  |2 

c  At rx  =  \urx  (ttx  +  At ix  j  +  vrxtp2  —  ulgt  (tfx  +  Atix  )| 

Define  ii,  (tlx  +  At ix)  =  Urx (0)  +  vja  +  vtgtAtix  - utgt (0)  - vtgrttx  - vtgtA tix 
=>  c2AT  =  (ur  (ttx  +  Atix)+ vtgtA trx ) •  (ur  (tlx  +  Atix )  +  vlgtAtrx ) 

=>  0  =  \nr  {ttx  +  Atix  )f  +  2“,  (ttx  +  Atix  )  •  a t„  +  (|v„  |2  -  c2 )  At)x 


Solving  the  quadratic  equation  for  tp  and  choosing  the  causal  solution  where  Atrx  >  0  : 
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«r(f*+A  tix>V„  + 

K*  [t,x )  = - 

This  is  the  general  case  where  the  transmitter  and  receiver  are  not  assumed  to  be  co¬ 
located  or  moving  with  the  same  speed. 

In  the  case  of  a  monostatic  radar,  ur  (ta  +  Atix )  =  —ut  (ttt  +  Atix )  and  vrx  =  vtx  .  This  yields 
Equation  (14)  from  Chapter  II. 

-u,  (ttx  +  A  tix  )-vtx  +  J(ur(ttx  +  Atix)-vtxf  +(c2  -|vte|2W(y|2 

(/tt )  =  7~  2\ 

(c  -H ) 

C.  DOPPLER-SHIFT  FOR  INCOMING  TARGET 

In  the  case  of  an  incoming  target, 

Ur  K ) '  \t  =  K  K  )| '  | \t  |  cosd  8°)  =  -  %  (ttx  )| '  \\, 

Starting  from  Equation  (13),  the  derivation  of  Equation  (22)  is  as  follows: 


_“t(0-Vtgt+^l 

\a,(Of  ■ 

t2  +(c2-'VIgt2)\ut(ttx)  I' 

(C2-  V 

i2 

*  ) 
tgt\  ' 

\llt(tJ\2 

z-\ut(t,x)  r  +c2\ut(ttx)\2 

SJ 

+ 

£ 

tsT 

l 

(c2- 

' \u,(t,x)\ 

vf) 

,  2 
(c  -yrgt 

) 

ri,(tlx)  (c-  v,.„ ) 
(c-|vw|Xc  +  |vv|) 
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(C  +  Kr) 


At  (t,  )  : 

nr'  tx' 


Equation  (23)  can  be  derived  from  Equations  (14)  and  (22)  similarly  as  follows: 

_  -  ft  (U  +  AU )  |  K I  +  i/ft  (U  +  Ar,t)|2  |vrv  |2  +  (c2  -  |v,r  |2 )  I U,  (ttx  +  Atix)\2 

,  2  I-  |2. 

(c  ~KI ) 

[Since  \vtx  |  =  0] 

[By  definition  of  ut(ttx)  and  since  IvJ  =  0] 


u,(ttx  +  At,x)\c 


~  A 4 


(ft(ob  V 


c  + 


tgt 


[Using  Equation  (17)  for  Atix] 


CM. 


■,(U)|  +  (|^|-|v|)ft(U)| 


(c+lvl) 


(c+  i 


tgt 


The  result  of  Equation  (24)  can  then  simply  be  shown  as  follows 


A'A)  =  Ar/A)+AUO  = 


2ft(u) 

(c  +  |uJ) 


A?/4  +  «AU)  = 


2ft(U-  +  »AU-) 

(c+ftpl) 


2(ft(U)|-«AUftgr|) 

(c+ftg.l) 


Atp(tix+nAttx)-Atp(ta)  = 


-2 

«AU 

(c  + 

) 
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APPENDIX  B.  MATLAB  CODES 


The  MATLAB  codes  used  in  for  this  study  are  listed  in  this  Appendix. 


A.  FUNCTIONS 

The  MATLAB  functions  written  for  the  model  are  listed  here. 

1.  Generation  of  Linear  Trajectory  Path 


function  [Pos]  =  linear_tra j ectory (P_start, Vel, time) 

%  This  function  is  to  calculate  the  position  of  an  object  at 
%  a  given  value  of  time  (or  a  row  vector  of  time  values) . 

%  Assumes  (for  the  moment)  that  the  object  is  travelling  linearly 
%  with  constant  speed,  Vel 

Q, 

O 

%  Inputs: 

%  -  Vel  =  speed  vector,  3x1  column  vector  of  [vx;vy;vz] 

%  -  P  start  =  Position  at  time  =  0,  3x1  column  vector  of  [x0;y0;z0] 
%  -  time  =  time  instants  to  compute  position  for  [any  row  vector] 

Q, 

O 

%  Outputs: 

%  -  P  =  Matrix  of  positions  [X;Y:Z]  corresponding  to  [time] 

Q, 

O 

%  By:  Tng  YS  (Jul  2012) 

[A1,B1]  =  size (time); 
if  (A1  ~=  1) 

error ('Only  single  values  or  row  vectors  allowed') 

end 

Pos  =  repmat ( P_start,  1 ,  B1 )  +  repmat (Vel, 1, Bl) . *repmat (time, 3,1); 


2.  Computation  of  Wave  Interception  Time 


function  [t_i]  = 

intercept_time (obj_start_pos, ob j_velocity, wave_start_pos ) 

o, 

o 

%  This  function  computes  the  time  that  a  wavefront  travelling  at 
%  speed  of  light  from  tx_pos  takes  to  intercept  an  object  tgt  starting 
from 

%  tx_pos. 

O, 

o 
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%  Inputs: 

%  -  obj  start_pos  =  Object  position  at  wave  transmission, 3xN  column 
vector 

%  -  obj_velocity  =  Obj  (linear)  velocity,  3xN  column  vector 
%  -  wave^start  pos  =  Wave  start  position,  3xN  column  vector 

O, 

o 

%  Outputs: 

%  -  t  i  =  lxN  row  vector  of  intercept  times 

O, 

o 

%  By:  Tng  YS  (Jul  2012) 

[A1,A2]  =  size (obj_start_pos) ; 

[A3,A4]  =  size (obj_velocity) ; 

[A5,A6]  =  size (wave_start_pos) ; 

if  (A1  ~=  3  |  A3  ~=  3  |  A5  ~=  3  | |  A4  ~=  1) 

error ('All  position  inputs  must  be  3xN  matrices,  velocity  is  3x1  ') 

end 

if  (A2  ~=  A6) 

error ( 'Dimensions  of  the  position  matrices  must  be  the  same.') 

end 

c  =  299792458;  %  m/s  -  EM  wave  speed  =  Speed  of  light  in  vaccuum 
u  =  ob j_start_pos  -  wave_start_pos ; 

%  This  assumes  constant  velocity.  If  variable  velocity,  change  function 
%  to  accept  the  velocity  vector  accordingly, 
v  =  repmat (obj_velocity,  1,  A6)  ; 

u_dot_v  =  dot(u,v); 
v_dot_v  =  dot(v,v); 
u_dot_u  =  dot(u,u); 
c2  minus_v2  =  cA2  -  v_dot_v; 

t  i  =  (u  dot  v  +  sqrt (  u  dot_v.A2  +  c2  minus  v2.*u  dot  u))./c2  minus  v2 
end 


3.  Generation  of  Rectangular  Pulse 


function  [t,f_t,mf]  =  rect_pulse ( time_vec, f c, f s , pw, PRI , t_start ) 

%  This  code  is  used  to  generate  the  values  of  a  single  (amplitude 
modulated) 

%  rectangular  pulse  transmit  waveform  at  the  time  instances  specified 
by  time  vector 

%  The  matched  filter  for  a  single  pulse  duration  is  also  returned 

Q, 

O 

%  The  waveform  is  described  by  a  carrier  frequency,  pulse  width, 
repetition 

%  period,  transmission  start  time  (relative  to  t  =  0) . 
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%  time  vec  is  the  vector  of  time  values  to  compute  waveform  for; 

%  fc  is  the  carrier  frequency 
%  pw  is  pulse  width 

%  PRI  is  the  interval  between  pulses 
%  t_start  is  the  start  transmision  time 

%  Calculate  length  of  t 
[Nrow,Ncol]  =  size (time  vec) ; 
num_samples  =  Ncol; 
if (Nrow  ==  1 ) 

num_samples  =  Ncol; 
elseif (Ncol  ==  1) 

num_samples  =  Nrow; 
time^vec  =  time  vec'  ; 

else 

error ('time  must  be  a  vector'); 

end 

check  time  vec (time  vec, t_start, pw, PRI ) ; 
t  =  double (time^vec) ; 
if  (pw  >=  PRI) 

%  This  becomes  a  CW  signal  . . . 
f  t  =  ones ( 1 , num_samples ) ; 
if  (t_start  >=  t(l)) 

start_ind  =  find(t  <  t_start) ; 
f_t (1 : start_ind)  =  0; 

end 

else 

%  Compute  the  intervals  where  the  pulses  occur  &  set  these  to  1 
f  t  =  zeros ( 1 , num_samples ) ; 
pulse_start  =  t_start; 
pulse_end  =  t_start  +  pw; 

start  ind  =  max (1, find (t  >=  pulse  start,  1,  '  first'  ))  ;  %  Find 
smallest  index 

if  (t(end)  >  pulse_end) 

end_ind  =  find(t  <=  pulse  end, 1 ,' last' )  ;  %  Find  largest  index 

else 

end  ind  =  num^samples; 

end 

f  t (start  indrend  ind)  =  1; 

end 

first  start  ind  =  start  ind; 

%  Generate  matched  filter  baseband  response 
mf  =  con j (f_t (first  start_ind : end  ind) ) ; 

Nmf_orig  =  length (mf);  Nmf  =  2 Anextpow2 (Nmf_orig) ; 
mf  =  [mf, zeros ( 1 , Nmf -Nmf  orig) ] ; 

%  Upconversion  to  RF 

f  t  =  f_t.*exp(li  *  2  *  pi  *  fc  *  (time_vec  +  t_start) ) ;  %  Multiply 
Baseband  with  Carrier 
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4, 


Generation  of  Barker- Coded  Binary  Phase  Modulated  Pulse 


function  [t,f_t,mf]  = 

BarkerCode (time_vec, fc, fs, codelength, subpw, PRI, t_start) 

%  This  code  is  used  to  generate  the  values  of  a  single  Barker  coded 
%  pulse  transmit  waveform  at  the  time  instances  specified  by  time 
vector 

%  The  matched  filter  for  a  single  pulse  duration  is  also  returned 

O, 

o 

%  The  waveform  is  described  by  a  carrier  frequency,  pulse  width, 
repetition 

%  period,  transmission  start  time  (relative  to  t  =  0) . 

%  time  vec  is  the  vector  of  time  values  to  compute  waveform  for; 

%  fc  is  the  carrier  frequency 
%  fs  is  the  sampling  frequency 

%  code  length  is  the  length  of  the  Barker  code  desired  (Up  to  13) 

%  subpw  is  sub-pulse  width 
%  PRI  is  the  interval  between  pulses 
%  t_start  is  the  start  transmision  time 

%  Calculate  length  of  t 
[Nrow,Ncol]  =  size (time  vec) ; 
num_samples  =  Ncol; 
if (Nrow  ==  1 ) 

num_samples  =  Ncol; 
elseif (Ncol  ==  1) 

num^samples  =  Nrow; 
time^vec  =  time  vec'  ; 

else 

error ('time  must  be  a  vector' ) ; 

end 

%  The  actual  pulse  width  is  7  sub-pulse  widths 
pw  =  codelength  *  subpw; 

check  time  vec (time  vec, t_start, pw, PRI ) ; 

t  =  double (time_vec)  ; 

switch (codelength) 
case  2 

Code  =  [+1,-1]; 
case  3 

Code  =  [+1, +1, -1] ; 
case  4 

Code  =  [+1, +1, -1, +1] ; 
case  5 

Code  =  [+1,+1,+1,— 1 , +1 , ] ; 

case  7 

Code  —  [+1,+1,+1,— 1,-1, +1,-1]; 

case  11 

Code  =  [+1 , +1 , +1 ,  — 1 , — 1 , — 1 , +1 , — 1 , -1 , +1 , — 1 ] ; 

case  13 
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end 


Code  =  [+1 , +1 , +1 , +1 , +1 , — 1 , — 1 , +1 , +1 , — 1 , +1 , — 1 , +1 ] ; 

otherwise 

error ( 'Not  a  valid  Barker  code' ) ; 


if  (pw  >=  PRI) 

error ('PW  cannot  be  greater  than  PRI  for  the  chosen  Barker  code') 

else 

f  t  =  zeros (l,num_samples) ; 

%  pulse_index_len  =  zeros ( 1 , codelength) ;  %  For  debug 
for  nn  =  1: codelength 

%  Compute  the  intervals  where  the  subpulses  occur  &  set  these 
to 

%  the  Barker  code 

pulse_start  =  t_start  +  (nn-1 ) *subpw; 
pulse  end  =  t_start  +  nn*subpw; 

start_ind  =  max (1, find (t  >=  pulse  start, 1 ,' first' )) ;  %  Find 
smallest  index 

if  (nn  ==  1) 

first_start  ind  =  start_ind; 

end 

if  (t(end)  >  pulse_end) 

end  ind  =  find(t  <=  pulse  end, 1 ,' last' ) ;  %  Find  largest 

index 

else 

end  ind  =  num_s ample s ; 

end 

%  pulse^ind  length  (nn)  =  end_ind  -  start__ind;  %  For  debug 
f_t ( start_ind : end_ind)  =  Code(nn); 

end 

end 

%  Generate  matched  filter  baseband  response 
mf  =  con j (f_t (first  start_ind : end  ind) ) ; 

Nmf  orig  =  length (mf);  Nmf  =  2 Anextpow2 (Nmf  orig)  ; 
mf  =  [mf, zeros (1, Nmf -Nmf  orig) ] ; 

%  Upconversion  to  RF 

f  t  =  f_t.*exp(li  *  2  *  pi  *  fc  *  (time^vec  +  t  start));  %  Multiply 
Baseband  with  Carrier 


5.  Generation  of  Linear-Frequency-Modulated  Pulse 

function  [t, f_t, mf , mf_freq_ind]  =  lfm(time  vec, . . . 

f c, f s , PulseBW, pw, PRI , t  start, direction) 

%  This  code  is  used  to  generate  the  values  of  a  single  Linear  FM 
%  rectangular  pulse  transmit  waveform  at  the  time  instances  specified 
by  time  vector 

%  The  matched  filter  for  a  single  pulse  duration  is  also  returned 

O, 

o 


85 


%  The  waveform  is  described  by  a  carrier  frequency,  pulse  width, 
repetition 

%  period,  transmission  start  time  (relative  to  t  =  0) . 


%  time  vec  is  the  vector  of  time  values  to  compute  waveform  for; 

%  fc  is  the  carrier  frequency 
%  PulseBW  is  the  chirp  BW 
%  pw  is  pulse  width 

%  PRI  is  the  interval  between  pulses 

%  direction  (  'up'  or  'down'  is  the  chirp  direction) 

%  t_start  is  the  start  transmision  time 

%  Set  chirp  direction 
switch  (direction) 
case  'down' 

downchirp  =  1; 
case  'up' 

downchirp  =  0; 
otherwise 

error (' Select  chirp  up  or  down'); 

end 

%  Calculate  length  of  t 
[Nrow,Ncol]  =  size (time  vec)  ; 
num^samples  =  Ncol; 
if (Nrow  ==  1 ) 

num_samples  =  Ncol; 
elseif (Ncol  ==  1) 

num_samples  =  Nrow; 
time^vec  =  time  vec' ; 

else 

error ('time  must  be  a  vector'); 

end 

check  time  vec (time  vec, t_start, pw, PRI ) ; 

t  =  double (time^vec)  ; 
mu  =  PulseBW/pw; 

if  (pw  >=  PRI) 

%  This  becomes  a  CW  signal  . . . 
f  t  =  ones ( 1 , num_samples ) ; 
if  (t_start  >=  t  ( 1 ) ) 

start_ind  =  find(t  <  t_start) ; 
f_t (1 : start_ind)  =  0; 

end 

else 

%  Compute  the  intervals  where  the  pulses  occur  &  set  these  to  1 
f  t  =  zeros (l,num_samples) ; 
pulse_start  =  t_start; 
pulse_end  =  t_start  +  pw; 

start  ind  =  max (1, find (t  >=  pulse  start,  1 ,' first' ))  ;  %  Find 
smallest  index 
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Find  largest  index 


if  (t(end)  >  pulse_end) 

end_ind  =  find(t  <=  pulse_end, 1 , ' last' ) ;  % 

else 

end  ind  =  num^samples; 

end 

%  Chirp  Up 

f  t ( start_ind : end_ind)  =  exp (li*pi*mu* (t (start  ind:end  ind) - 
t ( start_ind) ) .  A2 ) ; 

%  Chirp  Down 
if (downchirp) 

f  t (start  indrend  ind)  =  f  t(end  ind: -1: start  ind); 

end 

end 

%  Generate  the  unshifted  baseband  matched  filter  response  &  zero  pad  to 
%  length  that  is  power  of  2  (for  fft  &  easy  computation  of  freq  index) 
mf  =  con j ( f_t ( start_ind : end_ind) ) ; 

Nmf  orig  =  length (mf);  Nmf  =  2 Anextpow2 (Nmf^orig) ; 
mf  =  [mf, zeros (1, Nmf -Nmf  orig) ] ; 

%  This  is  the  rectangular  pulse  where  f  t  is  simply  the  amplitude 
modulation 

f  t  =  f_t.*exp(li  *  2  *  pi  *  fc  *  (time^vec  +  t  start));  %  Multiply 
Baseband  with  Carrier 


6.  Correlation  Receiver 

function  [corr  out  r , delay_ind, f req_ind]  ... 

corr  rx (w  r,m_^f,rng  width  desired, f_width  desired, fstep  desired, fs , Ts ) 

O, 

o 

%  w_r  is  the  received  waveform 
%  m  f  is  the  baseband  matched  filter 

%  f_width  desired  (Hz)  is  the  range  of  frequencies  to  generate  on  each 
side 

%  fstep  desired  (Hz)  is  the  desired  frequency  shift  between  each  row 
%  Ts  is  the  sampling  interval 


o, 

o 

'k-k-k-k'k-k-k-k'k-k-k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k-k-k'k-k'k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k-k-k'k-k-k-k'k'k 

%  Relevant  Physics  Constants  (from  NIST) 

g, 

o 

: k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

c  =  299792458;  %  m/s  -  EM  wave  speed  =  Speed  of  light  in  vaccuum 

Nmf  =  length (m  f ) ; 

Nwr  =  length (w  r)  ; 


o 

%  Generate  the  bank  of  matched  filters  matched  to  shifted  frequencies 
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%  Compute  number  of  range  cells  around  scene  center  to  compute  on  each 
side 

rng_sample  =  c*Ts/2;  %  Two-way  range  (Correct  for  monostatic) 

N  rng  cell  =  ceil (rng  width  desired/rng_sample) ;  %  cells 

fstep  =  fs/Nmf;%  Hz  per  frequency  index  step  in  FFT  output 
nshift  =  f step_desired/f step; 

Nf  =  ceil(f_width  desired/f step^desired) ;  %  No.  of  freq  rows  on  each 
side 

Nf  limit  =  5000;  %  To  limit  and  avoid  out  of  memory  problems 
if  (Nf  >  Nf_limit) 

Nf  =  Nf_limit; 

disp('***  WARNING:  Frequency  steps  truncated  to  5000  on  each  side 

*  *  *  '  )  • 
end 

mf  freq_vec  =  li*2*pi* [-Nf : 1 :Nf ] '  *nshift/Nmf ; 
index  vec  =  l-Nmf/2 :1 :Nmf/2; 
n  shift_mat  =  repmat (mf  freq_vec, 1 , Nmf ) ; 
index  mat  =  repmat (index  vec, 2*Nf+l, 1) ; 
f  shift_mat  =  exp(n  shift_mat . *index  mat); 
mf  mat  =  f  shift  mat . *repmat (m  f,2*Nf+l,l); 

mf_conj  =  conj (mf  mat);  %  MATLAB  dot  product  will  conjugate  it  again 


%  Generate  associated  frequency  indices 
freq_ind  =  [-Nf : 1 :Nf ] ' *nshif t*f step; 

%  Size  of  mf  mat  is  needed  later 
col  mf  =  Nmf;  row  mf  =  2*Nf+l; 

%  Ensure  that  received  signal  is  at  least  same  length  as  the  matched 
filter 

%  for  correlation  to  work  (This  can  happen  between  mf  is  always 
extended 

%  to  2AN  lengths  for  efficient  FFT) 
if (Nwr  <  col_mf ) 

w  r  =  [w_r,  zeros (l,col_mf-Nwr) ] ;  %  Append  zeros 

end 


'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Correlation  Receiver 

^  ■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k'k-k-k-k-k-k'k-k'k-k'k 

w  r2  =  repmat ([ zeros ( 1 , N  rng  cell),w  r, zeros (1,N  rng_cell ) ] , row  mf,l); 
[~,Nwr2]  =  size(w_r2); 

corr  out  r  =  zeros (row  mf,2*N  rng  cell); 
for  nn  =  1 : (2*N  rng  cell+1) 

mf2  =  [zeros (row__mf, (nn-1 )), mf_conj , zeros (row  mf , Nwr2-col_mf- (nn- 

1 )  )  ]  ; 

corr_out  r(:,nn)  =  dot(mf2,w  r2,2); 

end 

%  Generate  associated  time  delay  indices  (referenced  to  start  of  w_f) 
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Associated  delay  from 


delay  ind  =  [-N  rng  cell:l:N  rng  cell]*Ts;  % 
center 

%  Normalise  by  the  waveform  energy 

r_0  =  dot (m_f , m_f ) ; 

corr  out  r  =  corr^out_r/r_0; 

end 


7.  Create  the  Time  Map  for  Given  Target-Sensor  Geometry 

function  [time  map, img  map,x  ind, y_ind]  =  ... 

create_timemap ( scene_ctr , x_width, x_step, y_width, y__step, tx_pos , rx_pos ) ; 

%  Compute  a  reception  time  map  for  the  transmitter  &  reception  geometry 
%  based  on  cartesian  coordinates  around  a  scene  and  initialise  an  image 
map 

O, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Relevant  Physics  Constants  (from  NIST) 

Q. 

O 

-k  ' k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

c  =  299792458;  %  m/s  -  EM  wave  speed  =  Speed  of  light  in  vaccuum 

%  Form  N_r  x  N  xr  x  2  matrix  to  capture  x,y-coord  for  each  pt 
N  x  =  ceil (x_width/x_step) ; 

N_y  =  ceil (y_width/y_step) ; 
x  ind  =  [-N  x:N  x] *x_step+scene_ctr ( 1 ) ; 
y_ind  =  [-N_y:N_y] *y_step+scene_ctr (2) ; 
spatialmap  =  zeros (2*N_y+l,2*N  x+1,2); 
spatialmap ( : , : , 1 )  =  repmat (x  ind, 2*N_y+l, 1) ; 
spatialmap 2 )  =  repmat (y_ind' , 1,2 *N_x+l) ; 
img  map  =  zeros (2*N  y+l,2*N  x+1,1); 

%  For  each  point  on  spatial  map,  compute  the  propagation  time. 
tx2pos_x  =  spatialmap 1 )  -  tx_pos(l,:); 
tx2pos_y  =  spatialmap 2 )  -  tx_pos(2,:); 

tx2_pos  =  sqrt (tx2pos_x . A2  +  tx2pos_y . A2 ) ; 

rx2pos_x  =  spatialmap 1 )  -  rx_pos(l,:); 
rx2pos_y  =  spatialmap 2 )  -  rx_pos(2,:); 

rx2_pos  =  sqrt (rx2pos_x . A2  +  rx2pos_y . A2 ) ; 

time  map  =  (tx2_pos  +  rx2_pos)/c; 
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8. 


Update  the  Image  Map  based  on  Current  Range  Samples 


function  [img  map  updated]  =  update  img  map(img  map, time  map, f r , trx, Ts ) 

%  Updates  the  img  map  with  the  (a  single)  function  value  based  on 
%  given  reception  time  instant. 

temp  img  map  =  zeros (size (img^map) ) ; 

gl  =  find (time  map  <  (trx+Ts/2)); 
temp  img  map(gl)  =  fr; 
g2  =  find (time  map  <  (trx-Ts/2) ) ; 
temp  img  map(g2)  =  0; 

img  map  updated  =  img  map  +  temp  img  map; 


9.  Check  Time  Vector  for  Required  Length  Constraints  before  Processing 

function  check  time  vec(time  vec,t  start, pw, PRI ) ; 

%  Convert  back  to  double  precision  for  easy  processing 
%  **  Watch  out  for  lost  accuracy  beyond  16th  significant  digit  ! ! 

%  **  Use  only  the  actual  symbolic  time  vector  for  calculations  ! ! 
t  =  double (time^vec) ; 

%  Check  that  the  time  vector  contains  only  a  single  pulse 
if  (t(end)  >  t_start  +  PRI  I  I  t(end)  >  t(l)  +  PRI) 

error ('Time  vector  contains  more  than  1  PRI.  Please  separate  the 
pulses  for  calculations!'); 

%  Processing  the  PRI  separately  makes  matched  filtering  easy  to 
implement . 

%  This  is  what  is  done  in  practical  receiver  systems  as  well. 

end 

%  Check  that  time  vector  is  sufficient  to  complete  1  pulse 
if  (t(end)  <  t_start  +  pw) 

error ('Time  vector  contains  less  than  1  Pulse  Width.  Will  not  yield 
correct  matched  filtering!'); 
end 

%disp ( 'Checked  OK'); 
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B 


SCRIPTS 


This  Section  provide  examples  of  scripts  which  make  use  of  the  implemented  functions 
to  simulate  different  imaging  radar  configurations 

1.  Example  of  Multistatic  Radar 

a.  Scenario  Definition 

%  This  script  defines  the  target  characteristics 

%  It  should  be  called  for  all  M-files  to  ensure  a  consistent  sensor 
%  definition 

%  Simulation  parameters 

plot  scene  =1;  %  This  is  to  set  whether  to  do  the  quiver  plot  or  not. 


O, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Relevant  Physics  Constants  (from  NIST) 

O, 

o 

"k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

c  =  299792458;  %  m/s  -  EM  wave  speed  =  Speed  of  light  in  vaccuum 

o, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Imaging  Parameters 

Q, 

O 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Set  Region  of  Interest  around  pulse  and  scene 
scene_ctr  =  [-50,50,0]; 

rng  width  desired  =  50;  %  meters  (on  each  side) 
rng_step  =  0.1;  %  meters  (per  cell) 

rng  resol  desired  =10;  %  meters 

cr  rng  width  desired  =50;  %  meters  (on  each  side) 
cr_step  =0.1;  %  meters  (per  cell) 

cr_rng  resol_desired  =10;  %  meters  (on  each  side) [Not  used  in 
Multistatic] 

O, 

o 

i^^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Sensor  Waveform  &  Sampling  Time  Definition 

O, 

o 

"k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

fc  =  150e6;  %  Carrier  frequency  in  Hz 

fs  =  le9;  %  Sampling  frequency;  Should  be  at  least  Nyquist+ 

Ts  =  1/fs; 

PRF  =  0.1;  %  Hz 
PRI  =  1/PRF;  %  sec 
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wavef orm^choice  =  wavef orm^choice  vec(qq); 
switch (waveforra_choice) 

case  1  %  Rectangular  Pulse 
PW  =  le-7;  %  sec 
PulseBW  =  1/PW;  %  Hz 
case  2  %  Barker  Code 

SubPW  =  (2/3)*le-7;  %  sec 

codelength  =  13;  %  No.  of  sub-pulse  widths 
PulseBW  =  1/ SubPW;  %  Hz 
PW  =  codelength*SubPW;  %  sec 
case  3  %  Linear  FM 
PW  =  le-6;  %  sec 

PulseBW  =  c/ (2*rng_resol  desired) ;  %Hz 
lfm  direction  =  'up'  ; 

end 

%  For  matched  filter  generation 

rng_sample  =  c*Ts/2;  %  meters  (based  on  sampling  rate) 
f  width  desired  =  0*PulseBW;  %  Hz  (on  each  side  on  zero  Doppler  cut) 
f step_desired  =  ... 

max (10e3, round (f  width  desired/500 ) ) ;  %Hz  (difference  betw  each  row) 


o, 

o 

"k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Sensor  Position,  Velocity  Definition 

Q, 

O 

'k-k-k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k-k-k'k-k-k-k'k-k-k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k-k-k'k-k-k-k'k'k 

%  Notes: 

%  -  All  theta  are  measured  counter-clockwise  from  x-axis 
%  -  tx  range  is  distance  measured  from  sensor  to  scene  center 
%  -  velocity  bearing  =  0  =>  moving  away  from  scene  center  along  x 
%  -  velocity  bearing  =  180  =>  moving  away  from  scene  center  along  -x 

tx  pos  theta  =0;  %  deg  -  Angle  between  tx  &  scene  center  at  start 
tx  init  height  =0;  %  meters 
tx  init_range  =  100e3;  %  meters 

tx  init_pos  =  tx  init_range* [cosd (tx_pos  theta) ; . . . 

sind (tx_pos_theta) ; tx  init  height];  %meters  -  @  time  =  0; 
tx_spd  =0;  %  m/s 

tx  vel_theta  =  180;  %deg  -  velocity  bearing 

tx_vel  =  tx_spd* [cosd (tx_vel_theta) ; sind (tx_vel_theta) ; 0 ] ;  %  m/s 

%  Set  to  same  as  tx  for  1-channel  mono-static  radar.  Adjust  according 
%  otherwise 

rx  pos  theta  =  theta  vec(qq);  %  deg  -  Angle  between  tx  &  scene  center 
at  start 

rx  init  height  =0;  %  meters 

rx  init_range  =  rx  init  rng  vec(qq); 

rx  init_pos  =  rx  init^range* [cosd (rx_pos  theta) ; . . . 

sind (rx_pos_theta) ; rx_init  height];  %meters  -  @  time  =  0; 
rx_spd  =0;  %  m/s 

rx  vel_theta  =  180;  %deg  -  velocity  bearing 

rx^vel  =  rxspd* [cosd (rx_vel  theta) ; sind (rx_vel_theta) ; 0 ] ;  %  m/s 
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-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Target  Position,  Velocity  Definition 

Q, 

O 

%  Notes: 

%  -  All  theta  are  measured  counter-clockwise  from  x-axis 
%  -  velocity  bearing  =  0  =>  moving  away  from  scene  center  along  x 
%  -  velocity  bearing  =  180  =>  moving  away  from  scene  center  along  -x 

tgt_spd  =0;  %  m/s 

tgt_init_pos  =  [-50;50;0];  %meters  -  @  time  =  0; 

%savefile  id  =  [num2str(tgt  spd/le3) , ' k' ] ;  %  Used  in  script  for 
analysis 

tgt_vel_theta  =0;  %  deg  -  velocity  bearing 

tgt_vel  =  tgt_spd* [cosd (tgt_vel_theta) ; sind (tgt_vel_theta) ; 0 ] ;  %  m/s 


o, 

o 

%  Plot  the  relative  positions  and  vel  on  a  quiver  plot  if  selected 

o, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

if (plot_scene) 

scale  =  1000; 
v_scale  =  200; 

if (qq  ==  1) 

h  scene  =  figure; 

else 

figure (h^scene) ; 

end 

if  (qq  ==  1  &&  tgt_spd  ==  0) 

scatter (tgt_init_pos ( 1 ) /scale, tgt_init_pos (2) /scale,  60,  . . . 

'd' , ' Markerf acecolor '  ,  '  k'  )  ; 
elseif  (qq  ==  1) 

quiver (tgt_init_pos (1) /scale, tgt_init_pos (2) /scale, . . . 
tgt_vel ( 1 ) /v^scale, tgt_vel (2 ) /v_scale, ' MaxHeadSize' , 10 , . . . 
'LineWidth' , 3, ' color' , ' k' ) 
end 

hold  on 

if  (tx_spd  ==  0) 

scatter (tx_init_pos (1) /scale, tx_init_pos (2) /scale,  40,  .  .  . 

's',' Markerf acecolor '  ,  '  r '  )  ; 

else 

quiver (tx_init_pos (1) /scale, tx_init_pos (2) / scale, . . . 
tx_vel (1) /v^scale, tx  vel(2)/v  scale, ' MaxHeadSize' , 2 ,.. . 
'LineWidth' , 3, ' Color' , ' r' ) ; 

end 

if  (rx_spd  ==  0) 

scatter (rx_init_pos (1) /scale, rx_init_pos (2) /scale,  40,  .  .  . 

'o' , ' Markerf acecolor '  ,  '  r'  )  ; 

else 

quiver (rx_init_pos (1) /scale, rx_init_pos (2) /scale,  .  .  . 
tx  vel (1) /v^scale, tx  vel(2)/v  scale, ' MaxHeadSize' , 2 ,.. . 
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'LineWidth' , 3, ' Color' ,  ' r'  )  ; 

end 

hold  off; 

%  xlim ( [-5, 15] ) ;  ylim ( [-15, 15] ) ; 
grid  on 

xlabel ( 'X  (km) ' ) ; 
ylabel ( ' Y  ( km) ' ) ; 

end 


b.  Main  Script 

%  This  script  is  used  to  generate  the  multistatic  config  in  the  thesis 

o, 

o 

%close  all;  %  Close  all  figures  (These  take  up  memory  too!) 

format  long;  %  Set  MATLAB  display  format 

clock  =  tic;  %  To  keep  track  of  the  script's  run  time 

%  These  define  the  angles  and  ranges  for  the  receivers  (for  multistatic) 
theta_vec  =  [-45,0,30,90];  %  degrees 

rx_init_rng_vec  =  [sqrt (2) *100e3, 100e3, 100e3, 100e3*sind (30) ] ; 

wavef orm_choice_vec  =  [3, 3, 3, 3];  %  waveform  selection  for  each  tx/rx 

pair 

N  sys  =  length (theta  vec) ;  %  No.  of  systems  simulated  in  this  script 
for  qq  =  1 : N  sys 

scenario  definition  multistatic;  %  Call  scenario  file  (multistatic) 

if  (qq  ==  1) 

rng_sample  =  c*Ts/2; 

N  rng  cell  =  ceil (rng  width  desired/rng_sample) ;  %  cells 

range_prof ile  =  zeros (N  sys , 2 *N_rng_cell+l ) ; 

trx  range_prof ile  =  zeros (N_sys , 2 *N  rng_cell+l); 

range  alignment  ind  =  zeros (N  sys,l); 

tx_pos  =  zeros ( 3 , N_sys ) ; 

rx_pos  =  zeros ( 3 , N_sys ) ; 

end 


%  Create  the  time  range  to  simulate  based  on  PW 

^  'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k 

t  start  =  0;  %sec  -  Required  start  time  (integral  multiple  of  Ts  needed) 
t  end  =  t  start+PW+Ts;  %sec  -  Required  end  time ( integral  multiple  of  Ts 
needed) 

^  ^^icic^ic^-k-k-k'k-k-k-k'k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k'k-k-k-k'k-k-k-k'k-k-k-k-k-k-k-k'k-k-k-k'k-k-k-k-k-k-k-k'k 

%  Generation  of  Time  Sequences 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


n  =  round(t  start/Ts); 
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N_desired  =  round (t_end/Ts) ; 

%  For  actual  N,  adjust  number  of  samples  for  time  expansion  due  to 
motion 

N  =  ceil (c/ (c- (tx_spd+rx_spd+tgt_spd) ) *N_desired) ;  %  worst  case 
if  ( (N-n)  >  5e6+1000 ) 
n  samples  =  N-n 

fprintf ( [ 'Too  many  samples  simulated.  Memory  problems  may  occur. 

r  .  .  . 

'It  is  recommended  to  reduce  sequence  size  to  <  5  million  ', . 
'samples .  \n']); 

fprintf ([' For  longer  sequences,  alternative  is  to  simulate  as  '... 

'separate  vectors  and  combine  subsequently  if  needed.  \n' ] ) ; 

A1  =  0; 

while  (A1  ~=  1  &&  A1  ~=  2) 

A1  =  input ( 'Continue  ?  (1  to  continue,  2  to  stop)  :  '); 

if  (A1  ==  2) 

return  %  exit  the  script 

end 

end 

end 

%  Create  evenly  spaced  (@  Ts)  reference  time  vector 

%  --  Time  vector  starts  from  nTs  and  ends  at  N*Ts  to  allow  flexibility 
%  in  the  start  transmission  time 
%  --  Length  is  N-n+1  samples 
ttx  =  n*Ts : Ts : N*Ts ;  %  sec 

%  Compute  propagation  time 

[trx, tp, t_ctr]  =  total_propagation_time (tgt_init_pos, tgt_vel, . . . 
tx  init_pos,tx  vel,rx  init_pos,rx  vel,ttx); 

S ^  -k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Recreate  evenly  spaced  receive  samples  using  search 

-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Desired  outcome 
trx^desired  =  tp(l)  +  ttx; 
trx  original  =  trx; 
ttx  original  =  ttx; 
tp_original  =  tp; 
t_ctr_original  =  t_ctr; 

%  Search  algorithm 

alpha  =  0.99;  %  This  controls  how  fast  search  converges 
for  curr  sample  =  2:N+1 
curr^sample 

%  Initialise  before  entering  while  loop 
desired  time  =  trx  desired (curr^sample) ; 
actual_time  =  trx (curr^sample) ; 
diff  from  desired  =  actual_time-desired  time; 

ttx (curr  sample)  =  ttx (curr  sample)  -  alpha  *  diff  from  desired; 


iterations  =  0;  %  This  is  to  prevent  infinite  loops 
threshold  =  max ( le-1 6*desired  time, . . . 

le-18);%  le-16  is  the  limit  of  double  precision 

while  ( (abs (diff_from  desired)  >  threshold)  &&  (iterations  <  1000)) 

iterations  =  iterations  +  1; 

%  Compute  propagation  time 

[trx (curr^sample) , tp (curr_sample) , t_ctr (curr_sample) ] . . . 

total_propagation_time (tgt_init_pos, tgt_vel, tx_init_pos, . . . 

tx  vel,rx  init_pos,rx  vel, ttx (curr  sample)); 

%  Compute  difference  from  desired 

actual  time  =  trx (curr  sample); 

dif f_f rom^desired  =  actual__time-desired_time; 

%  Update  for  next  iteration 

ttx (curr  sample)  =  ttx (curr  sample)  -  alpha  *  diff  from  desired 

end 

end 


^  "k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Create  Transmitted  Waveform  &  Matched  Filter 

^  'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

switch (waveform_choice) 

case  1  %  Rectangular  Pulse 

[~,f_t,mf]  =  rect_pulse (ttx_original- 

ttx^original (1) , fc, fs, PW, . . . 

PRI , t_start) ; 
case  2  %  Barker  Code 

[~,f  t,mf]  =  BarkerCode (ttx^original-ttx_original (1) , fc, fs, . . . 
codelength, SubPW, PRI , t_start)  ; 
case  3  %  Linear  FM 

[~,f  t,mf]  =  lfm(ttx  original-ttx  original (1) , fc, fs, .. . 
PulseBW, PW, PRI , t  start, lfm_direction) ; 

end 

S ^  'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Propagation 

switch (waveform_choice) 

case  1  %  Rectangular  Pulse 

[~,f_r,~]  =  rect_pulse (ttx-ttx (1) , fc, fs, PW, PRI, 0) ; 
case  2  %  Barker  Code 

[~,f_r,~]  =  BarkerCode (ttx-ttx (1) , fc, fs, .. . 
codelength, SubPW, PRI, 0) ; 
case  3  %  Linear  FM 

[~,f_r,~]  =  lfm (ttx-ttx ( 1 ), fc,  fs ,..  . 

PulseBW, PW, PRI , 0 , lfm  direction) ; 


end 
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^  -k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Reception 

9>  ■k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Demodulate  the  carrier  frequency  using  the  original  time  vector 
mixer_phase  =  0;  %  Assumes  Ideal  mixer  (Perfect  phase  &  timing) 
f  r  demod  =  f  r . *exp (-li*2*pi*fc*ttx) *exp (-li*mixer_phase) ; 

%  Sampling/Decimation/Downsampling  (if  required) 

sample  start  =  1;  %  To  allow  sampling  at  different  start  samples 
N  dec  =1;  %  This  is  the  number  of  times  to  downsample, 
m  f  =  mf ( sample_start : N  dec : end); 
f  r  demod  =  f^r^demod ( sample  start :N  dec : end); 

Ts  =  Ts*N_dec; 
fs  =  fs/N_dec; 

%  Correlation  receiver 

[corr  out  r,  delay  ind,mf  freq_ind]  .  .  . 

=  corr  rx (f_r_demod,mf , rng_width  desired,  f  width  desired,  .  .  . 
f step_desired, fs, Ts)  ; 

%  Select  range  profile  (Max  is  selected) 

[Nrow,Ncol]  =  size  (corr  out^r) ; 
zero^cut  =  floor (Nrow/2+1 ) ; 

range_prof ile (qq, : )  =  corr_out_r ( zero_cut, : ) ;  %  Zero-Doppler  cut 
range  alignment  ind(qq)  =  find_range_alignment(trx-t_ctr(l),Ts); 

trx_range_prof ile (qq,  :  )  =  delay_ind  +  trx(l); 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Image  Reconstruction 

^  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

%  Create  Time  Map 

tx_pos ( : , qq)  =  tx_init_pos  +  t_start*tx_vel; 
rx_pos  (  :  ,  qq)  =  rx_init_pos  +  t_start*rx__vel ; 

[time  map, img  map,x  ind, y_ind]  =  ... 

create  timemap ( scene_ctr , rng  width  desired,  rng  step,  .  .  . 

cr_rng_width_desired, cr_step, tx_pos ( : , qq) , rx_pos ( : , qq) ) ; 

%  Check  and  Update  Time  Map  for  each  sample  in  range  profile 
Max  trx  =  max (max (time  map));%  Max  trx  in  desired  scene 
Min_trx  =  min (min (time  map));%  Min  trx  in  desired  scene 
for  11  =  1 : Ncol 

if ( (trx_range_prof ile (qq, 11)  >  Min  trx)  &&  ... 

(trx  range_profile (qq, 11)  <  Max  trx)) 

[img  map]  =  update  img  map (img  map, time  map, . . . 

range_prof ile (qq, 11) , trx_range_prof ile (qq, 11) , Ts) ; 

end 

end 

if (qq  ==  1) 


97 


img  map^cum  =  zeros (size (img  map)); 

end 

img  map^cum  =  img  map_cum+img_map; 
eval(['img  map  ' , num2str (qq)  ,  '  =  img  map; ' ] ) ; 
end 

h  img2  =  surf  (x  ind,  y__ind,  abs  ( img  map  cum) /N  sys, ' edgecolor' , ' none' ) ; 
xlabel ( 'X  (meters) ' ) ; 
ylabel ( 'Y  (meters) ' ) ; 

elasped_time  =  toe (clock) 


2.  Example  of  SAR 

a.  Scenario  Definition 


This  script  defines  the  target  characteristics 

It  should  be  called  for  all  M-files  to  ensure  a  consistent  sensor 
definition 


%  Simulation  parameters 
plot  scene  =  1;  %  This  is 


to  set  whether  to  do  the  quiver  plot  or  not. 


O, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Relevant  Physics  Constants  (from  NIST) 

O, 

o 

'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

c  =  299792458;  %  m/s  -  EM  wave  speed  =  Speed  of  light  in  vaccuum 


O, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


%  Imaging  Parameters 

g, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  For  moving  SAR 
SAR_speed  =  50.0;  %m/s 
SAR  halfwidth  =  1.5e3;  %  meters 
SAR^downrange  =  10e3;  %  meters 

%  Set  Region  of  Interest  around  pulse  and  scene 
scene_ctr  =  [0,200,0]; 

rng  width  desired  =  100;  %  meters  (on  each  side) 
rng_step  =  2;  %  meters  (per  cell) 

cr_rng_width  desired  =  100;  % 


meters  (on  each  side) 
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cr^step  =  2;  %  meters  (per  cell) 


o, 

o 

'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

%  Sensor  Waveform  &  Sampling  Time  Definition 

g, 

o 

-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

fc  =  100e6;  %  Carrier  frequency  in  Hz 

fs  =  1000e6;  %  Sampling  frequency;  Should  be  at  least  Nyquist+ 

Ts  =  1/fs; 

wavelength  =  c/fc;  %  m 
rng_resol  desired  =10;  %  meters 

cr_rng  resol  =  wavelength/ . . . 

(2*(atan(SAR  halfwidth/SAR  downrange) ) ) ;  %  meters  (on  each  side) 

waveform^choice  =  3; 
switch (waveform_choice) 

case  1  %  Rectangular  Pulse 
PW  =  le-7;  %  sec 
PulseBW  =  1/PW;  %  Hz 
case  2  %  Barker  Code 

SubPW  =  2/3*le-8;  %  sec 

codelength  =  13;  %  No.  of  sub-pulse  widths 
PulseBW  =  1 / SubPW;  %  Hz 
PW  =  codelength*SubPW;  %  sec 
case  3  %  Linear  FM 
PW  =  le-6;  %  sec 

PulseBW  =  c/ (2*rng  resol  desired) ;  %Hz 
lfm  direction  =  'up' ; 

end 

Nyquist_sampling  =  wavelength/2; 

PRI  =  Nyquist  sampling/SAR  speed;  %  sec 
PRF  =  1 /PRI ;  %  Hz 

%  For  matched  filter  generation 

rng_sample  =  c*Ts/2;  %  meters  (based  on  sampling  rate) 
f  width_desired  =  0*PulseBW;  %  Hz  (on  each  side  on  zero  Doppler  cut) 
f step_desired  =  ... 

max (10e3, round (f  width  desired/500));  %Hz  (difference  betw  each  row) 


g, 

o 
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%  Sensor  Position,  Velocity  Definition 

g, 

o 
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%  Notes: 

%  -  All  theta  are  measured  counter-clockwise  from  x-axis 
%  -  tx  range  is  distance  measured  from  sensor  to  scene  center 
%  -  velocity  bearing  =  0  =>  moving  away  from  scene  center  along  x 
%  -  velocity  bearing  =  180  =>  moving  away  from  scene  center  along  -x 

tx  init  height  =0;  %  meters 


99 


tx  init_pos  =  [SAR  downrange; -SAR  halfwidth; . . . 

tx  init  height];  %meters  -  @  time  =  0; 

tx_spd  =  SAR_speed; 

tx_vel  = [ 0 ; SAR_speed; 0 ] ;  %  m/s 

%  Set  to  same  as  tx  for  1-channel  mono-static  radar.  Adjust  according 
%  otherwise 

rx  init  height  =  tx  init  height;  %  meters 
rx  init_pos  =  tx  init_pos; 
rx_spd  =  tx_spd; 
rx  vel  =  tx  vel;  %  m/s 


O, 

o 
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%  Target  Position,  Velocity  Definition 

O, 

o 

%  Notes: 

%  -  All  theta  are  measured  counter-clockwise  from  x-axis 
%  -  velocity  bearing  =  0  =>  moving  away  from  scene  center  along  x 
%  -  velocity  bearing  =  180  =>  moving  away  from  scene  center  along  -x 

tgt_spd  =0;  %  m/s 

tgt_init_pos  =  [0;0;0] ;  %meters  -  @  time  =  0; 
tgt_vel_theta  =0;  %  deg  -  velocity  bearing 

tgt_vel  =  tgt_spd* [cosd (tgt_vel_theta) ; sind (tgt_vel_theta) ; 0] ;  %  m/s 


%  Plot  the  relative  positions  and  vel  on  a  quiver  plot  if  selected 


if (plot_scene  ==  1) 
scale  =  1000; 
v_scale  =  60; 
h  scene  =  figure; 
if  (tgt_spd  ==  0) 

scatter (tgt_init_pos ( 1 ) /scale, tgt_init_pos (2) /scale, 60, . . . 
'd' , ' Markerf acecolor '  ,  '  k'  ) ; 

else 

quiver (tgt_init_pos (1) /scale, tgt_init_pos (2) /scale, . . . 


tgt_vel ( 1 ) /v_scale, tgt  vel (2 ) /v^scale, ' MaxHeadSize' , 15 , ' LineWidth' , 3 ) 
end 

hold  on 

if  (tx_spd  ==  0) 

scatter (tx_init_pos (1) /scale, tx_init_pos (2) /scale, 40, . . . 

'd' , ' Markerf acecolor '  ,  '  k'  ) ; 

else 

quiver (tx_init_pos ( 1 ) /scale, tx_init_pos (2) /scale,  .  .  . 
tx  vel ( 1 ) / v^scale,  tx  vel (2)/v  scale,'  MaxHeadSize'  ,  15 ,  .  .  . 
'LineWidth' , 3, ' Color' , ' r' ) ; 

end 
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hold  off; 

xlim ( [-1, 11] ) ;  ylim ( [-2,2] ) ; 
grid  on 

xlabel ( 'X  (km) ' ) ; 
ylabel ( ' Y  ( km) ' ) ; 

end 


b.  Main  Script 

%  This  script  is  used  to  simulate  a  moving  SAR  config  for  the  thesis 

O, 

o 

%close  all;  %  Close  all  figures  (These  take  up  memory  too!) 

format  long;  %  Set  MATLAB  display  format 

clock  =  tic;  %  To  keep  track  of  the  script's  run  time 

plot  scene  =0;  %  This  is  to  set  whether  to  do  the  quiver  plot  or  not. 

scenario_def inition_SAR;  %  Call  scenario  definition  file 

plot  scene  =0;  %  This  is  to  set  whether  to  do  the  quiver  plot  or  not. 

N_pulses  =  ceil (2*SAR  halfwidth/wavelength*2 ) ; 
t_start_vec  =  (0 : (N_pulses-1) ) *PRI; 
for  qq  =  l:N_pulses 

qq 

scenario  definition  SAR;  %  Call  scenario  definition  file 
if  (qq  ==  1) 

rng_sample  =  c*Ts/2;  %  Two-way  range  (Correct  for  monostatic) 

N  rng  cell  =  ceil (rng  width  desired/rng^sample) ;  %  cells 

range_prof ile  =  zeros (N_pulses, 2*N_rng_cell+l ) ; 

range  alignment  ind  =  zeros (N_pulses , 1 ) ; 

tx_pos  =  zeros ( 3 , N_pulses ) ; 

rx_pos  =  zeros ( 3 , N_pulses ) ; 

trx  range_prof ile  =  zeros (N_pulses, 2*N  rng  cell+1); 

end 


%  Create  the  time  range  to  simulate  based  on  PW 
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t_start  =  t_start_vec (qq) ;  %sec  -  Required  start  time  (integral 
multiple  of  Ts  needed) 

t  end  =  t  start+PW+Ts;  %sec  -  Required  end  time ( integral  multiple  of  Ts 
needed) 
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%  Generation  of  Time  Sequences 
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n  =  round (t_start/Ts) ; 

N_desired  =  round (t_end/Ts) ; 

%  For  actual  N,  adjust  number  of  samples  for  time  expansion  due  to 
motion 

N  =  ceil (c/ (c- (tx_spd+rx_spdttgt_spd) ) * (N_desired-n) ) ;  % 
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worst  case 


if  ( (N-n)  >  5e6+1000 ) 
n  samples  =  N-n 

fprintf ( [ 'Too  many  samples  simulated.  Memory  problems  may  occur. 

A 

f  •  •  • 

'It  is  recommended  to  reduce  sequence  size  to  <  5  million  ', . . . 
'samples .  \n' ] ) ; 

fprintf ([' For  longer  sequences,  alternative  is  to  simulate  as  '... 

'separate  vectors  and  combine  subsequently  if  needed.  \n' ] ) ; 

A1  =  0; 

while  (A1  ~=  1  &&  A1  ~=  2) 

A1  =  input ( 'Continue  ?  (1  to  continue,  2  to  stop)  :  '); 

if  ( A1  ==  2) 

return  %  exit  the  script 

end 

end 

end 

%  Create  evenly  spaced  (@  Ts)  reference  time  vector 

%  --  Time  vector  starts  from  nTs  and  ends  at  N*Ts  to  allow  flexibility 
%  in  the  start  transmission  time 
%  --  Length  is  N-n+1  samples 
ttx  =  n*Ts : Ts : (N+n) *Ts;  %  sec 

%  Compute  propagation  time 

[trx, tp, t_ctr]  =  total_propagation_time (tgt_init_pos, tgt_vel, . . . 
tx  init_pos,tx  vel , rx_init_pos , rx  vel,ttx); 
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%  Recreate  evenly  spaced  receive  samples  using  search 
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%  Desired  outcome 
trx_desired  =  tp(l)  +  ttx; 
trx  original  =  trx; 
ttx_original  =  ttx; 
tp_original  =  tp; 
t_ctr  original  =  t  ctr; 

%  Search  algorithm 

alpha  =  0.99;  %  This  controls  how  fast  search  converges 

for  curr  sample  =  2: (N+l) 

%  Initialise  before  entering  while  loop 
desired__time  =  trx  desired  (curr_sample)  ; 
actual_time  =  trx (curr  sample); 
diff_from  desired  =  actual_time-desired_time; 

ttx (curr  sample)  =  ttx (curr  sample)  -  alpha  *  diff  from  desired; 

iterations  =  0;  %  This  is  to  prevent  infinite  loops 
threshold  =  max (le-16*desired  time, . . . 

le-18);%  le-16  is  the  limit  of  double  precision 

while  ( (abs (diff  from  desired)  >  threshold)  &&  (iterations  <  1000)) 
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iterations  =  iterations  +  1; 

%  Compute  propagation  time 

[trx (curr^sample) , tp (curr^sample) , t_ctr (curr_sample) ] . . . 

=  total_propagation_time (tgt_init_pos, tgt_vel, . . . 
tx  init_pos , tx^vel , rx  init_pos,rx  vel, ttx (curr  sample)); 

%  Compute  difference  from  desired 

actual  time  =  trx (curr  sample); 

dif f_f rom_desired  =  actual_time-desired_time; 

%  Update  for  next  iteration 

ttx (curr  sample)  =  ttx (curr  sample)  -  alpha  *  diff  from  desired 

end 

end 
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%  Create  Transmitted  Waveform  &  Matched  Filter 
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switch (waveform_choice) 

case  1  %  Rectangular  Pulse 

[~,f_t,mf]  =  rect_pulse (ttx_original, fc, fs, PW,  .  .  . 

PRI , t_start) ; 
case  2  %  Barker  Code 

[~,f_t,mf]  =  BarkerCode (ttx_original, fc, fs,  .  .  . 
codelength, SubPW, PRI , t_start) ; 
case  3  %  Linear  FM 

[~,f  t,mf]  =  lfm(ttx  original, fc, fs, . . . 

PulseBW, PW, PRI , t  start, lfm  direction); 

end 
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%  Propagation 
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switch (waveform_choice) 

case  1  %  Rectangular  Pulse 

[~,f_r,~]  =  rect_pulse (ttx-ttx (1) , fc, fs, PW, PRI, 0) ; 
case  2  %  Barker  Code 

[~,f_r,~]  =  BarkerCode (ttx-ttx (1) , fc,  fs,  ..  . 
codelength, SubPW, PRI, 0) ; 
case  3  %  Linear  FM 

[~,f_r,~]  =  lfm (ttx-ttx ( 1 ), fc, fs ,.. . 

PulseBW, PW, PRI , 0 , lfm  direction) ; 

end 
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%  Reception 
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%  Demodulate  the  carrier  frequency  using  the  original  time  vector 
mixer_phase  =  0;  %  Assumes  Ideal  mixer  (Perfect  phase  &  timing) 
f  r  demod  =  f_r . *exp (-li*2*pi*fc*ttx) *exp (-li*mixer_phase) ; 
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%  Correlation  receiver 

[corr  out  r , delay_ind, mf_f req_ind]  . . . 

=  corr  rx (f_r_demod, mf , rng_width  desired, f  width  desired, . . . 
f step_desired, fs,  Ts)  ; 

%  Select  range  profile 
[Nrow,Ncol]  =  size (corr  out  r) ; 
zero  cut  =  f loor (Nrow/2+1 ) ; 

trx_range_prof ile (qq,  :)  =  delay^ind  +  trx(l)  -  t_start; 

range  alignment  ind(qq)  =  find  range  alignment (trx-t_ctr (1) , Ts) ; 

%  Based  on  Soumekh  (pp  359) 

%  Reverse  baseband  conversion  -  Note  this  requires  sampling  rate  to 
%  be  much  higher  that  the  carrier  frequency 
range_prof ile (qq, : )  =  corr^out  r(zero  cut,:).*... 
exp (li*2*pi*fc* (delay_ind+t_start)  )  ; 
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%  Image  Reconstruction 
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%  Create  Time  Map 

tx_pos ( : , qq)  =  tx_init_pos  +  t_start*tx_vel; 
rx_pos ( : , qq)  =  rx  init_pos  +  t  start*rx  vel; 

[time  map, img  map,x  ind, y_ind]  =  ... 

create  timemap ( scene_ctr , rng  width  desired,  rng  step,  .  .  . 

cr_rng_width_desired, cr_step, tx_pos ( : , qq) , rx_pos ( : , qq) )  ; 

%  Check  and  Update  Time  Map  for  each  sample  in  range  profile 
Max  trx  =  max (max (time  map));  %  Max  trx  in  desired  scene 
Min_trx  =  min (min (time  map));  %  Min  trx  in  desired  scene 
for  11=1: Ncol 

if ( (trx  range_profile (qq, 11)  >  Min  trx)  &&  ... 

(trx_range_profile (qq, 11)  <  Max^trx) ) 

[img  map]  =  update  img  map (img  map, time  map, . . . 

range_prof ile (qq, 11) , trx_range_prof ile (qq, 11) , Ts) ; 

end 

end 

if (qq  ==  1) 

[Nrow  img, Ncol  img]  =  size (img_map) ; 
img  map^cum  =  zeros (Nrow  img, Ncol  img) ; 

end 

img  map^cum  =  img  map_cum+img  map; 
end 
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%  Analysis  and  Graph  plotting 
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o 


h  img  =  surf (x  ind, y_ind, abs ( img  map  cum) /N_pulses, ' edgecolor' ,' none' ) 
xlabel ( 'Down  Range  (meters)');  ylabel ( 'Cross  Range  (meters)'); 
view ([0,90]); 
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if  (plot_scene  ==  1) 
figure (h  scene) 
hold  on 

quiver ( (tx_init_pos (1) +tx_vel (1) *trx (end) ) /scale,  .  .  . 
(tx_init_pos (2 ) +tx_vel (2 ) *trx (end) ) /scale, . . . 
tx_vel ( 1 ) /v_scale, tx_vel (2 ) /v_scale,  .  .  . 
'MaxHeadSize' , 2, ' LineWidth'  ,  3,  '  Color'  ,  '  r'  )  ; 
hold  off 

end 

elasped_time  =  toe (clock) 


3.  Script  for  Re-imaging  the  Scene  After  Generation 


%  This  script  to  allow  processing  the  range  profile  matrix  with  a 
different 

%  scene  of  interest  following  completion  of  the  original  run 

[Nrow,Ncol]  =  size (range_prof ile) ; 

%  Set  Region  of  Interest  around  pulse  and  scene 
scene_ctr  =  [0,200,0]; 

rng_width  desired  =50;  %  meters  (on  each  side) 
rng_step  =1;  %  meters  (per  cell) 

cr  rng  width  desired  =  600;  %  meters  (on  each  side) 
cr_step  =  1;  %  meters  (per  cell) 

%  For  movie  creation 

mov(l:Nrow)  =  struct ( 'edata' ,  [ ] , ' colormap' ,  []); 

for  qq  =  1 : Nrow 

qq 
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%  Image  Reconstruction 
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%  Create  Time  Map 

[time  map, img  map,x  ind, y_ind]  =  ... 

create  timemap ( scene  ctr,rng  width  desired, rng  step, . . . 

cr_rng_width_desired, cr_step, tx_pos ( : , qq) , rx_pos ( : , qq) ) ; 

%  Check  and  Update  Time  Map  for  each  sample  in  range  profile 
Max  trx  =  max (max (time  map));  %  Max  trx  in  desired  scene 
Min_trx  =  min (min (time  map));  %  Min  trx  in  desired  scene 
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for  11  =  1 : Ncol 

if ( (trx_range_prof ile (qq, 11)  >  Min  trx)  &&  ... 

(trx  range_prof ile (qq, 11)  <  Max^trx) ) 

[img  map]  =  update  img  map(img  map, time  map, . . . 

range_prof ile (qq, 11) , trx_range_prof ile (qq, 11)  ,  Ts)  ; 

end 

end 

if (qq  ==  1) 

img  map_cum  =  zeros (size (img  map) ) ; 

end 

img  map^cum  =  img  map_cum+img  map; 

h  img2  =  surf  (x  ind,  y__ind,  abs  ( img  map  cum)  /Nrow,  '  edgecolor'  ,  '  none'  ) 

xlabel ( 'X  (meters) ' ) ;ylabel ( 'Y  (meters) ' ) ; 

mov(qq)  =  getframe (gcf ) ;  %  For  movie  creation 

pause ( . 03 )  ; 

end 

%  For  movie  creation 

movie2avi (mov,  'SAR50  tgtlO  90deg#l . avi' ,  'compression',  'None' ) ; 

%  Display  the  image  map 

h  img3  =  surf (x  ind, y_ind, abs ( img  map  cum) /Nrow, ' edgecolor' ,' none' ) 
xlabel ( 'Down  Range  (meters)');  ylabel ( 'Cross  Range  (meters)'); 
view ([0,90]); 
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APPENDIX  C.  STUDY  ON  FLOATING  POINT  PRECISION 


This  appendix  is  mainly  to  highlight,  for  reference  purposes,  the  effect  of 
computational  precision  on  MATLAB’s  simulation  outputs  when  comparing  the 
Symbolic  Toolbox’s  Variable  Precision  Arithmetic  (VPA)  computation  (25  digits)  versus 
double  precision  floating  point  calculations. 

All  plots  presented  in  this  Appendix  show  the  trend  of  the  differences  between  the 
sample  interval  at  transmission  versus  the  interval  between  the  corresponding  samples,  as 
noted  at  reception,  for  a  simulation  scenario.  As  explained  in  Chapter  n,  the  spread  in  the 
samples  is  due  to  the  spread  in  propagation  time  due  to  transmitter  and  target  motions. 

In  these  plots,  a  linear  line  formed  by  extrapolating  from  the  first  two  samples  is 
subtracted  from  actual  difference  (as  compared  to  the  transmitted  sampling  time)  at 
reception  to  accentuate  any  non-linearity  present  in  the  trend  (i.e.,  if  the  spread  is  linear), 
then  the  deviation  from  a  linear  line  should  be  consistently  zero ,  or  else  the  differences 
from  the  linear  line  will  show  up  clearly. 

A  good  example  of  a  sanity  check  would  be  a  case  where  a  clear  1-D  Doppler 
shift  effect  can  be  predicted:  The  target  starts  from  [0,0,0]  km  with  velocity  [0,50,0]  m/s 
and  the  transmitter  starts  from  [0,50,0]  km  with  velocity  [0,-200,0]; 


Transmission  time  (ps) 


107 


This  gives  a  first  hint  that  the  limits  of  double  precision  have  been  reached,  since  the 
VPA  result  is  around  10-40  (nearly  zero)  but  the  double  precision  results  give  around  10- 

17 


In  the  next  example,  the  target  starts  from  [0,0,0]  km  with  velocity  [50,0,0]  m/s 
and  the  transmitter  starts  from  [500,0,0]  km  with  velocity  [100,-200,0].  A  sample  size  of 
1000  samples  is  always  assumed,  with  sampling  frequency  equal  to  2fc  where  the  value  of 
fc  is  given  in  each  graph.  Her e/c  acts  as  a  scaling  factor  for  the  deviation  from  a  linear 
line,  since  the  amount  of  deviation  from  a  linear  line  increases  with  decreased  fc  (i.e., 
increased  the  time  between  samples). 


Deviation  from  linear  for  lOps  pulse  (fc  =  0.1MHz) 


-  Symbolic  VPA 
 Double-Precision 


Transmission  time  (ps) 


Transmission  time  (ps) 


The  results  clearly  indicate  that  as  frequency  increases  (i.e.,  interval  between 
samples  shrink),  the  double  precision  floating  point  calculations  can  no  longer  accurately 
present  the  results  (due  to  insufficient  dynamic  range),  whereas  the  VPA  calculations 
continue  to  show  the  same  trend.  At  fc  greater  than  1  MHz  in  this  example,  it  can  be  seen 
that  the  floating  point  calculations  are  unable  to  accurate  depict  the  non-linearity  at  all, 
whereas  calculations  using  VPA  continue  to  faithfully  reproduce  the  non-linear  curve 
correctly. 
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x  10 


Deviation  from  linear  for  lOps  pulse  (fc  =  1MHz) 


Transmission  time  (ps) 


Deviation  from  linear  for  lOps  pulse  (fc  =  10MHz) 


While  accuracy  is  improved,  the  tradeoff  is  that  VPA  calculations  take  up  much 
more  computation  time,  hence  judicious  design  of  the  MATLAB  model  has  to  be 
undertaken  to  optimise  simulation  time  if  VPA  is  needed.  In  this  thesis  study,  the  use  of 
VPA  is  avoided  since  the  intent  is  to  study  fast  moving  targets  and  longer  pulses,  both  of 
which  do  not  require  the  precision  of  VPA  to  represent  accurately. 
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